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Abstract 


A new n um erical framework for solving conservation laws is being developed. This new 
approach differs substantially in both concept and methodology from the well-established 
methods— i.e., finite difference, finite volume, finite element, and spectral methods. It is 
conceptually simple and designed to overcome several key limitations of the above tradi- 
tional methods. 

A two-level scheme for solving the convection-diffusion equation 

du/dt + cl du/dx — fi d 2 u/dx 2 = 0 

is constructed and used to illuminate major differences between the current method and 
those mentioned above. This explicit scheme, referred to as the a-p scheme, has the unusual 
property that its stability is limited only by the CFL condition, i.e., it is independent of 
p. Also it will be shown that the amplification factors of the a-p scheme are identical 
to those of the Leapfrog scheme if p = 0, and to those of the DuFort-Frankel scheme if 
a = 0. These coincidences are unexpected because the a-p scheme and the above classical 
schemes are derived from completely different perspectives, and the a-p scheme does not 
reduce to the above classical schemes in the limiting cases. 

The a-p scheme is extended to solve the 1— D time-dependent Navier-Stokes equations 
of a perfect gas. Stability of this explicit solver also is limited only by the CFL condition. 
In spite of the fact that it does not use (i) any techniques related to the high-resolution 
upwind methods, and (ii) any ad hoc parameter, the current Navier-Stokes solver is capable 
of generating highly accurate shock tube solutions. Particularly, shock discontinuites can 
be resolved within one mesh interval. 

The inviscid (/i = 0) a-fi scheme is neutrally stable, i.e., free from numerical diffusion. 
Such a scheme generally can not be extended to solve the Euler equations. Thus, the 
inviscid version is modified. Stability of this modified scheme, referred to as the a-e scheme, 
is limited by the CFL condition and 0 < e < 1 where e is a special parameter that controls 
numerical diffusion. Moreover, if e = 0, the amplification factors of the a-e scheme are 
identical to those of the Leapfrog scheme, which has no numerical diffusion. On the other 
hand, if e = 1, these amplification factors unexpectedly become identical to each other and 
to the amplification factor of the highly diffusive Lax scheme. Note that, because the Lax 
scheme is very diffusive and it uses a mesh that is staggered in time, a two- level scheme 
using such a mesh is often associated with a highly diffusive scheme. The a-e scheme, 
which also uses a mesh staggered in time, demonstrates that it can also be a scheme with 
no numerical diffusion. 

The a-e scheme is extended to become an Euler solver. The extension has stability 
conditions similar to those of the a-e scheme. It also has the unusual property that numer- 
ical diffusion at all mesh points can be controlled by a set of local parameters. Moreover, 
it will be shown that the Euler extension is capable of generating accurate shock tube 
solutions with the CFL number ranging from 0.88 to 0.022. 
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1. Introduction 


The method of space-time conservation element and solution element [1-3] is a new 
numerical framework for solving conservation laws. This new approach differs substantially 
in both concept and methodology from the well-established methods— i.e., finite difference, 
finite volume, finite element, and spectral methods [4-8]. It is conceived and designed to 
overcome several key limitations of the above traditional methods. Thus, we shall begin this 
paper with a discussion of several considerations that motivate the current development. 

(a) A set of physical conservation laws is a collection of statements of flux conservation in 
space-time. Mathematically, these laws are represented by a set of integral equations. 
The differential form of these laws is obtained from the integral form with the assump- 
tion that the physical solution is smooth. For a physical solution in a region of rapid 
change (e.g., a boundary layer), this smoothness assumption is difficult to realize by 
a numerical approximation that can use only a limited number of discrete variables. 
This difficulty becomes even worse in the presence of discontinuities (e.g., shocks). 
Thus, a method designed to obtain numerical solutions to the differential form with- 
out enforcing flux conservation is at a fundamental disadvantage in modeling physical 
phenomena with high-gradient regions. Particularly, it may not be used to solve flow 
problems involving shocks. Contrarily, a numerical solution obtained from a method 
that also enforces flux-conservation locally (i.e., down to a computational cell) and 
globally (i.e.. over the entire computational domain) will always retain the basic phys- 
ical reality of flux conservation even in a region involving discontinuities. For this 
reason, the enforcement of both local and global flux conservation in space and time 
is a tenet in the current development. As will be shown, the concept of conservation 
element is introduced to serve this purpose. 

Among the traditional methods, finite difference, finite element, and spectral 
methods are designed to solve the differrential form of the conservation laws. Note 
that the set of integral equations usually solved in a finite-element scheme is equiv- 
alent to the differential form of the conservation laws assuming certain smoothness 
conditions. However, these integral equations generally are different from the integral 
equations representing the conservation laws. Even if they are cast into a conservative 
form, the resulting flux-conservation conditions generally do not represent the physical 
conservation laws. 

The finite volume method is the only traditional method designed to enforce flux 
conservation. A finite- volume scheme may enforce flux conservation in space only, or 
in both space and time. As a preliminary to this enforcement, a flux must be assigned 
at any interface separating two neighboring conservation cells. In a typical finite- 
volume scheme, it is evaluated by extrapolating or interpolating the mesh values at 
the neighboring cells. This evaluation generally requires an ad hoc choice of a special 
flux model among many models available [9-11]. Generally numerical results obtained 
are dependent on which model one chooses. 

Contrarily, by design, making the above ad hoc choice is not needed in the current 
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numerical framework. As will be shown, by using the concept of solution element , 
flux evaluation at an interface becomes an integral part of the solution procedure and 
requires no interpolation or extrapolation. 

(b) The numerical variables used in a spectral method, i.e., the expansion coefficients, 
are global parameters pertaining to the entire computational domain. As a result, a 
spectral method generally (i) lacks local flexibility and thus may be applied only to 
problems with simple geometry, and (ii) is hindered by the fact that it must deal with 
a full matrix that is diffucult to invert. 

By design, only local parameters will be used in the current method. Moreover, 
the set of discrete variables in any one of the numerical equations to be solved is either 
associated with a single solution element or a few immediately neighboring solution 
elements. Thus, one needs only to deal with a very sparse matrix. As will be shown, 
the maximum number of solution elements involved in a numerical equation of the 
current framework is independent of the order of accuracy of a particular scheme. 
Contrarily, the order of accuracy of a classical finite- difference scheme generally can 
be increased only by using variables of more mesh points in each of its equations. 
Usually, a side effect of this practice is an increase in numerical diffusion, a subject 
to be discussed shortly. Note that, in the absence of body force, direct physical 
interactions occur only among the immediate neighbors. The current design is also 
consistent with this physical reality. 

(c) Space and time traditionally are treated separately in the time marching schemes. 
Generally one obtains a system of ordinary differential equations with time being the 
independent variable after a spatial discretization. As an example, elements in the 
finite element method usually are used for spatial discretization. These elements are 
domains in space only. 

Because flux conservation is fundamentally a property in space- time, space and 
time are unified and treated on the same footing in the current method. For example, 
conservation elements and solution elements used in the time- dependent version of the 
current method are domains in space-time. The significance of this unified approach 
cannot be overemphasized. As will be shown, it makes it easier for a numerical 
analogue to share the same space-time symmetry of the physical laws. 

(d) In a finite-difference scheme, derivatives at mesh points are expressed in terms of mesh 
values of dependent variables by using finite-difference approximations. Accuracy of 
these approximations, especially those of higher-order accuracy, generally is excellent 
as long as dependent variables vary slowly across a mesh interval. However, it may 
not be adequate if these variables vary too rapidly. Thus, in a high-gradient region, 
e.g., a boundary layer, accuracy may demand the use of an extremely fine mesh. In 
turn, a prohibitively high computing cost may result. 

The current method avoids the above pitfall by expressing the numerical solution 
within a solution element as an expansion in terms of certain base functions. As in a 
spectral method, the expansion coefficients are considered as the independent numer- 
ical variables to be solved simultaneously . For simplicity, Taylor’s expansions will be 
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used in the current paper. For this special case, the expansion coefficients are inter 
preted as the numerical analogues of the derivatives. Note that (i) van Leer [12] also 
has attempted to improve accuracy by introducing two independent numerical vari- 
ables for each independent physical variable, and (ii) the current solution procedure 
has no resemblance with those used in compact difference schemes. 

(e) With a few exceptions, numerical diffusion generally appears in a numerical solution 
of a time-marching problem. In other words, the numerical solution dissipates faster 
than the corresponding physical solution. For a nearly inviscid problem, e.g., flow with 
a high Reynolds number, this could be very serious because numerical dissipation may 
overwhelm physical dissipation and cause a complete distortion of solutions. One may 
argue that numerical diffusion can be reduced by increasing the order of accuracy of 
the scheme used. However, because the order of accuracy of a scheme is generally 
determined with the aid of Taylor’s expansion, and the latter is valid only for a 
smooth solution, it has meaning only for a smooth solution. Thus the use of a scheme 
of higher-order accuracy may not reduce numerical diffusion associated with high- 
frequency Fourier components of a numerical solution. This is the reason that the 
Leapfrog scheme, which is free from numerical diffusion, can outperform schemes with 
higher-order accuracy in solving some wave equations [13]. 

In a study of finite-difference analogues of a simple convection equation [14], it was 
shown that a numerical analogue will be free from numerical diffusion if it does not 
violate certain space-time invariant properties of the convection equation. In other 
words, numerical diffusion may be considered as a result of symmetry- breaking by 
the numerical scheme. Because of its intrinsic nature of space-time unity, the current 
framework is an excellent vehicle for constructing a numerical analogue that shares 
the same space-time invariant properties with the physical equation. 

It is recognized that a certain amount of numerical diffusion may be needed 
to prevent large dispersive errors [15] that are often caused by the presence of high- 
frequency disturbances (such as round-off errors). Therefore, in the current paper we 
shall construct a model scheme for a simple convection equation in which its numerical 
diffusion is controlled by a single adjustable parameter. The numerical diffusion is shut 
off when this parameter is set to zero. Furthermore, an Euler solver will be constructed 
such that its numerical diffusion at ail mesh points can be controlled by a set of local 
parameters. 

(f) High- resolution upwind methods for solving the Euler equations [8], which we consider 
to be a branch of the finite volume method, are heavily dependent on characteristics- 
based techniques. For the 1-D time- dependent case, the characteristics are curves in 
space-time, and the coefficient matrix associated with the Euler equations [16] also 
can be diagonalized easily. As a result, these techniques are easy to apply. However, 
for multi-dimensional cases, the characteristics are 2-D or 3-D surfaces in space-time 
[17]. Moreover, the coefficient matrices cannot be diagonalized simultaneously by the 
same matrix [16]. Because of the above complexites, application of these techniques 
to multi-dimensional problems is much more diffucult. Furthermore, high-resolution 
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methods generally require the use of ad hoc parameters, e.g., flux-limiters and/or 
slope-limiters, and other ad hoc techniques. These ad hoc techniques generally are 
also difficult to extend to a space of higher dimension. 

Because the current framework is developed to solve multi-dimensional prob- 
lems, simplicity and generality weigh heavily in its design. Thus, we do not use 
characteristics- based techniques, and also try to avoid using ad hoc techniques. More- 
over, the concept of characteristics generally is not applicable to the Navier-Stokes 
equations, which are non-hyperbolic in nature. Therefore, the above decision also 
makes it easier for the current framework to solve the Navier-Stokes equations. 

This completes the discussion of the motivation for the current development. In sum- 
mary, the development is guided by the following requirements: (i) To enforce both local 
and global flux consevation in space and time with flux evaluation at an interface being 
an integral part of the solution procedure and requiring no interpolation or extrapolation; 
(ii) To use local discrete variables such that the set of variables in any one of the numer- 
ical equations to be solved is associated with a set of immediately neighboring cells; (iii) 
Space and time are unified and treated on the same footing; (iv) Mesh values of depen- 
dent variables and their derivatives are considered as independent variables to be solved 
simultaneously; (v) To minimize numerical diffusion, a numerical analogue should be con- 
structed, as much as possible, to be compatible with the space-time invariant properties 
of the corresponding physical equations; and (vi) To exclude the use of the characteristics- 
based techniques, and to avoid the use of ad hoc techniques as much as possible. It is 
the purpose of this paper to show that the above requirements can be met with a simple 
unified numerical framework. 

For any reader who is interested in getting an advance idea on how simple the present 
method can be, he is referred to the computer program listed at the end of the present 
paper. It is a shock-tube-problem solver constructed using the present method. The 
simplicity of the solver is easily appreciated by a comparison of the listed program and a 
typical program associated with high-resolution upwind methods. Not only is the listed 
program much smaller in size (it is self-contained and the main loop contains only 33 lin es), 
but it contains no Fortran statements such as “if”, “amax”, and “a min” which are used so 
often in the programs implementing high-resolution methods. The absence of the above 
Fortran statements in the listed program results from the effort in avoiding the use of the 
ad hoc techniques in the development of the present method. In spite of its simplicity, it 
will be shown in Sec. S that the present solver is capable of generating highly accurate 
shock tube solutions. 
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2. The a-fi Scheme 


In this section, we consider a dimensionless form of the 1-D convection diffusion equa- 


tion, i.e., 


du du d 2 u 

~§i + a di-' i d^- 


( 2 . 1 ) 


where the convection speed a, and the viscosity coefficient p (> 0) are constants. Let 
xi = x, and x 2 = t be considered as the coordinates of a two-dimensional Euclidean 
space E 2 . By using Gauss’ divergence theorem in the space-time E 2 , it can be shown that 
Eq. (2.1) is the differential form of the integral conservation law 



• ds = 0. 


( 2 . 2 ) 


As depicted in Fig. 1, here (i) S(V ) is the boundary of an arbitrary space-time region V 
in E 2 , (ii) h = (au- pdu/dx, u ) is a current density vector in E 2 , and (iii) ds = dan with 
da and n, respectively, being the area and the outward unit normal of a surface element 
on S(V). Note that (i) h ■ ds is the space-time flux of h leaving the region V through the 
surface element ds, and (ii) all mathematical operations can be carried out as though E 2 
were an ordinary two-dimensional Euclidean space. 

At this juncture, note that the conservation law given in Eq. (2.2) is formulated in 
a form in which space and time are unified and treated on the same footing. This unity 
of space and time is also a tenet in the following numerical development. It is a key 
characteristic that distinguishes the current method from most of the traditional methods. 

Let ft denote the set of mesh points (j,n) in E 2 (dots in Fig. 2(a)) where n = 
0, ±1/2, ±1, ±3/2, ±2, ±5/2, . . ., and, for each n, j = n ± 1/2, n ±3/2, n ±5/2, . . .. There is 
a solution element (SE) associated with each (;, n) G ft. Let the solution element SE(j, n) 
be the interior of the space-time region bounded by a dashed curve depicted in Fig. 2(b). It 
includes a horizontal line segment, a vertical line segment, and their immediate neighbor- 
hood. For the following discussions, the exact size of this neighborhood does not matter. 

For any (x,t) G SE(j,n), u(x,t), and h(x,t), respectively, are approximated by 
u*(x,t;j,n) and h*(x,t;j,n) which we shall define shortly. Let 

«•(*,* ;j,n) = u] + (u x )](x - Xj) + (U t ) 7(t - n (2.3) 

where (i) u], (u x )], and (u t )] are constants in SE(j,n), and (ii) (. Xj,t n ) are the coordinates 
of the mesh point ( j,n ). Note that 


u *( x j, t n ',j, n) = u", 


^u*(i,t; j,n) 

dx 


= (Ux)”, 


du*(x,t;j,n) 

dt 


= (■ ut )”• (2-4) 


Moreover, if we identify u”, (u z )” , and (u t )", respectively, with the values of u . dujdx , 
and du/dt at (xj,t n ), the expression on the right side of Eq. (2.3) becomes the first -order 
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Taylor’s expansion of u(x, t) at (xj,t n ). As a result of these considerations, uj, (u r )", and 
(u t )” will be considered as the numerical analogues of the values of u, du/dx, and du/dt 
at ( Xj,t n ), respectively. 

We shall require that u = u*(x,t ;j,n) satisfy Eq. (2.1) within SE(j,n). As a result 
of Eq. (2.4), this implies that 

( u t)j = ~ a {u x )j. ( 2 . 5 ) 

Because Eq. (2.3) is a first-order Taylor’s expansion, the diffusion term in Eq. (2.1) has no 
counterpart in Eq. (2.5). As a result, the diffusion term has no impact on how u*(x, t ; j, n ) 
varies with time within SE(j, n). However, as will be shown shortly, through its role in the 
numerical analogue of Eq. (2.2), it does influence time-dependence of numerical solutions. 
Note that, for a higher-order scheme, how u*(x,t ;j,n) varies with time within SE (j,n) 
will be influenced by the presence of the diffusion term. Combining Eqs. (2.3) and (2.5), 
one has 


u*(x,t;j,n) = u" +(u r )? [(r-Xj)-a(f-f n )], (x, t) G SE(j, n). (2.6) 

Because h = ( au — fidujdx^u)^ we define 

h*(x,t;j,n) = (au*(x,t;j,n) - ndu*(x,t-,j,n)/dx, u*(x,t;j,n)). (2.7) 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 2(a)) referred to 
as conservation elements (CE’s). As depicted in Figs. 2(c) and 2(d), a CE with its top- 
right (top-left) vertex being the mesh point (j, n) G 0 is denoted by CE_(j, n) (CE+(j, n)). 
Obviously the boundary of CE (CE+(j,n)), excluding two isolated points B and C 
(C and D), is formed by the subsets of SE (j, n) and SE (j - 1/2, n - 1 /2) (SE(j + 1/2 , n - 
1/2)). The current approximation of Eq. (2.2) is 

F±(j, n) d = <f h* ■ ds = 0 (2.8) 

JS(CE ± U,n)) 

for all (/, n) G fi. In other words, the total flux leaving the boundary of any conservation 
element is zero. Note that the flux at any interface separating two neighboring CE’s 
is calculated using the information from a single SE. As an example, the interface AC 
depicted in Figs. 2(c) and 2(d) is a subset of SE(/, n). Thus the flux at this interface is 
calculated using the information associated with SE (j,n). Also note that an SE is the 
interior of a space-time region. _Thus the vertices B, C, and D, strictly speaking, do not 
belong to any SE. As a result, h * is not defined at these points. However, contributions 
to the above integral from these isolated points are zero no matter what values of h* are 
assigned to them. For this reason, one may simply exclude them from the above surface 
integration. 

Because the surface integration across any interface separating two neighboring CE’s is 
evaluated using the information from a single SE, obviously the local conservation condition 
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Eq. (2.8) will lead to a global conservation relation, i.e., the total flux leaving the boundary 
of any space-time region that is the union of any combination of CE s will also vanish. 

Because each S(CE±(j, n)) is a simple closed curve in E 2 (see Fig. 1), the surface 
integration in Eq. (2.8) can be converted into a line integration. Let 

g* d = (— u*, au* - pdu*/dx), and dr d = (dx,dt) (2.9) 

Thus, dr is normal to ds and points in the tangential direction of the line segment joining 
the two points (x, t ) and (x + dx,t + dt). Because ds = ±(cft, —dx) [1, p.14], we have 

h* -ds = ±g* -dr (2.10) 


where the upper (lower) sign should be chosen if the 90° rotation from ds to dr is in the 
counterclockwise (clockwise) direction. By combining Eqs. (2.8) and (2.10), one concludes 
that 


F±(j,n) = j>s 


dr 


S(C£±0,n)) 


( 2 . 11 ) 


Note that the notation c.c. indicates that the line integration should be carried out in the 
counterclockwise direction. Substituting Eq. (2.6) into Eq. (2.11), and using the fact that 
the boundary of a CE is formed by the subsets of two SE’s, one has 


(ax) 2 


F±(j, n) = ±(1/2) (1 — v 2 + + (1 ~ 1/2 — 0( u x)j± 1/2 


+ 


2(1 =F 1 /) 


AX 


(<- 


u 


-l/2\ 
j±l/2 ) 


( 2 . 12 ) 


where 


v 


def CLAt 

— ? 
AX 


and 



4/iAt 

(ax) 2 


(2.13) 


Note that (i) the parameter v is the Courant number, and (ii) a more efficient method of 
flux evaluation will be presented later in this section. 

With the aid of Eqs. (2.8) and (2.12), u] and (u*)” can be solved in terms of 

and ( u x)’j±ij 2 if 1 — 1/2 ± ^ ^ 0, i.e., for all SE(j, n), 


^ (j, n) = Q+ q(j — 1/2, n — 1/2) + Q- q(j + 1/2, n — 1/2) (1 - v 2 + £ ^ 0). (2.14) 


Here 


— / • \ def 

q(j,n) = 




(ax/4)(u x )" 


(2.15) 
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for all ( j , n) G ft, and 


<?+ = (1/2) 


and 


Q- = (1/2) 


1 + u 

1 - i/ 2 + ( 

1 - V 
1 -I/ 2 


1 — v 2 — £ 

-(1 - iQ(l - ^ - Q 

i - 1/ 2 + e 

-(1 + t/)(l - Z/ 2 - 0 

1 - I' 2 + f 


(2.16) 


(2.17) 


1 - v 2 + £ 

Because numerical variables at a higher time level can be evaluated in terms of those at a 
lower time level by using Eq. (2.14), it defines a marching scheme. Furthermore, because 
this scheme models Eq. (2.1) which is characterized by two parameters a and /i, hereafter 
it will be referred to as the a- fj. scheme. 

As a preliminary for future developments, we apply Eq. (2.14) successively and obtain 


q(j, n + 1) = ( Q+) 2 q(j - 1, n) + ( Q+Q _ + Q-Q+)q(j, n)+(Q-) 2 q (j + 1, n) 

( l -^ 2 + ^ 0 ). 


(2.18) 


A result of Eq. (2.18) is 


q(j, n + 1) -+ q(j,n) as At -+ 0, (2.19) 

if a, n, and ax are held constant. The proof follows from the fact that 

(<?+) 2 — ►0, (Q+Q- + Q-Q+) -* 1 and (Q_) 2 — ► 0, as At -»• 0, (2.20) 

if a, /z, and ax axe held constant. 

Alternatively, Eq. (2.19) can be proved using the fact that the total flux of h* leaving 
the boundary of any space-time region that is the union of any combination of CE’s 
vanishes. Consider the union of CE + (j, n + 1) and CE_(j + 1/2, n + 1 /2) (see Fig. 2). This 
union is a rectangle with the vertices (j + 1/2, n + 1), (j,n + 1), (;',n), and (j + 1/2, n). 
The flux leaving this rectangle through its two vertical edges approaches zero as At — ► 
0. Because the total flux leaving its boundary vanishes, one concludes that the total 
flux leaving its two horizontal edges also approaches zero as At — > 0. In other words, 
the flux entering the rectangle through the lower horizontal edge approaches that leaving 
through the upper horizontal edge as At — > 0. Because these two fluxes are evaluated using 
q(j, n) and q ( j , n+ 1), respectively, the above limiting condition implies a limiting relation 
between q(j,n ) and q(j, n + 1). Similarly, by considering the union of CE ~(j,n + 1) and 
CE+ (j — 1/2, n + 1/2), one obtains another limiting relation for q(j,n ) and q(j, n -f 1). 
Eq. (2.19) is a result of the above two limiting relations. 
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The a-fi scheme has several nontraditional features. They are summarized in the 

following remarks: 

(a) Space and time are unified and treated on the same footing in the construction of the 
a-/i scheme. 

(b) The expansion coefficients it” and (it r )” in Eq. (2.6) are treated as independent vari- 
ables, i.e., (u x )j is not expressed in terms of u"’s by using a finite-difference approxi- 
mation. 

(c) As a result of Eq. (2.12), each of the conservation conditions F±(j, n) — 0 involves 
only numerical variables associated with two neighboring SE’s. This fact remains true 
for a scheme of higher-order accuracy in which Eq. (2.3) is replaced by a Taylor’s 
expansion of higher-order. The contrast with the finite difference method and its 
physical significance were dicussed in Sec. 1. 

(d) The a-fi scheme has the simplest stencil, i.e., a triangle with a vertex at the upper time 
level and the other two vertices at the lower time level. Eq. (2.14), which relates nu- 
merical variables at these vertices, was derived using the flux conservation conditions 
F±(j , n) = 0. Because the flux at an interface separating two neighboring CE’s is eval- 
uated using information of a single SE, no interpolation or extrapolation is required. 
Moreover, accuracy of flux evaluation is enhanced by requiring that u — u* (x, t ; n) 
satisfy Eq. (2.1) within SE(j,n). This makes the use of characteristics-based tech- 
niques less necessary. 

(e) The a-/x scheme uses a mesh that is staggered in time. As will be explained in Ap- 
pendix A, for a two- level scheme using such a mesh, e.g., the Lax scheme [4, p.97], 
generally the numerical variable at (j, n + 1) does not approach that at (j, n) as At — ► 0, 
if a, n, and ax are held constant. This is a key reason why the Lax scheme is very dif- 
fusive when the Courant number v is small. According to Eq. (2.19), the a-/i scheme 
is an exception to the above general rule. 

(f) Eq. (2.1) can be solved numerically using the Leapfrog/DuFort-Frankel scheme [4, 
p.161]. This scheme is reduced to the Leapfrog scheme [4, p.100] if diffusion is absent 
(i.e., n = 0), and to the DuFort-Frankel scheme [4, p.114] if convection is absent (i.e., 
a = 0). It is well known that a solution of any of the above schemes is formed by 
two decoupled solutions with each being associated with a mesh that is also staggered 
in time. Traditionally the von Neumann stability analysis for the above schemes is 
performed without taking into account this decoupled nature [4]. In Appendix A, it 
is performed for each decoupled solution using the mesh depicted in Fig. 2(a). It is 
shown that the amplification factors of the Leapfrog/DuFort-Frankel scheme are 


A± 


£ cos(#/2) - iv sin(0/2) ± y^cos(0/2) - *Vsin(0/2)] 2 + 1 - f 2 

___ 


( 2 . 21 ) 


Here the amplification factors are defined to be those between the time levels n and 
n + 1, i.e., they axe the amplification factors of the solution after two marching steps. 
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The reason behind this definition is that the mesh points at the time levels n and 
n + 1 axe not staggered. Hereafter the same definition will be used for other schemes. 
Let 1 — v 2 ^ 0. Then the amplification factors of the current a-p scheme (see 
Eq. (6.9)) are identical to those given by Eq. (2.21) except that the parameter £ should 

be replaced by f = f £/(l — v 2 ). Because (i) £ = £ = 0 if p = 0, and (ii) v = 0 and 

<■» /i \ 

thus f = f, if a = 0, one concludes that are completely identical to those of 
the Leapfrog scheme if p = 0, and to those of the DuFort-Frankel scheme if a = 0. 
These coincidences are unexpected because the a-p scheme and the above classical 
schemes axe derived from completely different perspectives. Moreover, the a-p scheme 
is a two- level scheme with two variables u” and (u x )" associated with the mesh point 
(j, n), while the above classical schemes axe three-level schemes with a single variable 
u” associated with the same point. 

Because the amplification factors of the inviscid a-p scheme (i.e., the a-p scheme 
with p — 0) axe identical to those of the Leapfrog scheme, the former, as in the case 
of the latter, is neutrally stable (i.e., free of numerical diffusion) if u 2 < 1. Note that 
the case with /i = 0 and i/ 2 = 1 is ruled out by the assumption 1 — i/ 2 -f £ ^ 0 of 
Eq. (2.14). Similarly, the pure- diffusion a-p scheme (i.e., the a-p scheme with a = 0), 
as in the case of the DuFort-Frankel scheme, is unconditionally stable. Furthermore, 
it is proved in Sec. 6 that the stability of the general a-p scheme, as in the case of the 
Leapfrog/DuFort-Frankel scheme, is independent of /x, and restricted only by the CFL 
condition, i.e., u 2 < 1. The a-/x scheme is the only two-level explicit scheme known 
to the author to possesss the above properties. Also it will be shown later that the 
same stability condition is retained by a natural 1-D time- dependent Navier-Stokes 
extension of the a-p scheme. 

Because stability of the a-p scheme is restricted only by the CFL condition, the 
stability bound for At is proportional to ax. In contrast, the stability condition of 
a typical classical explicit scheme generally is more restrictive than the CFL condi- 
tion. For a small mesh Reynolds number, the stability bound for At is approximately 
proportional to (ax ) 2 for the MacCormack scheme [4, p.102]. 

Because a neutrally stable numerical analogue of the pure convection equation 


du du 

a +<, & =0 


( 2 . 22 ) 


usually becomes unstable when it is applied to a nonlinear inviscid generalization of 
Eq. (2.22), the inviscid a-p scheme will be modified in Sec. 3 such that it can be 
extended to model the Euler equations. In this new version, numerical diffusion is 
introduced in a way that allows its magnitude to be adjusted by a special parameter, 
(g) The conservation relations for CE + (/ — 1/2, n + 1/2) and CE_(j + 1/2, n + 1/2) (see 
Figs. 2(e) and 2(f)) are 


F+(j - 1/2, n + 1/2) = 0, and F-(j + 1/2, n + 1/2) = 0, (2.23) 
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il 


respectively. Combining Eqs. (2.12) and (2.23), and assuming 1 - v 2 - 0, one has 

q(j, n ) = Q+ q(j + 1/2, n + 1/2) + Q- q(j - 1/2, n + 1/2) (1 - u 2 - ^ 0 ). (2.24) 


Here 


and 


/ 


Q+ = ( 1 / 2 ) 


l + v 
1-u 2 



-(1 — v 2 + Z) \ 

-(1 -*/)(!- I/ 2 + Q ’ 
1 - v 2 - £ / 


Q- = (1/2) 


l-i/ 

-(1 ~ ^ 2 ) 
i - 1/ 2 - e 


l-^ + C \ 

-(1 + »/)(! - t ' 2 + o • 

i — i/ 2 — ^ / 


(2.25) 


(2.26) 


Eq. (2.24) defines a backward marching scheme, i.e., the numerical variables at the 
time level n axe determined in terms of those at the time level (n -f 1/2). Recall that 
both the forward marching scheme Eq. (2.14) and the backward marching scheme 
Eq. (2.24) are derived using the same set of conservation relations. As a matter of 
fact, Eqs. (2.14) and (2.24) are equvalent if (1 - v 2 ) 2 / (£) 2 is assumed. For the 
above reason, the a-y scheme may be referred to as a two-way marching scheme. For 
the case y > 0, it will be proved in Sec. 6 that the a-y scheme cannot be stable for 
both the forward and the backward marching directions except for the singular case 
v 2 -\ which is also on the threshold of instability. Thus, for all practical purposes 
the viscous a-y scheme is irreversible in time. On the other hand, it is neutrally stable 
for both the forward and backward marching directions, and thus is reversible in time, 
if y = 0, and v 2 < 1. Again, the a-y scheme is the only two-level explicit two-way 
marching scheme known to the author. 

(h) Several invariant properties of Eq. (2.1) with respect to space and time are discussed in 
[14]. In the same paper, these properties are also defined for the numerical analogues 
of Eq. (2.1). It is also shown that the neutral stability of several finite-difference 
analogues of Eq. (2.22) can be established by using their invariant properties with 
respect to space-time inversion. Because solutions of Eq. (2.22) do not dissipate with 
time, it is not surprising that solutions of a numerical analogue also will not dissipate 
with time, i.e., the scheme is neutrally stable, if it shares with Eq. (2.22) some space- 
time invariant properties. It will be shown in a future paper that the a-y scheme 
shares with Eq. (2.1) the same space-time invariant properties. Also note that these 
invariant properties are closely linked with the other properties discussed in (a), (e), 
(f), and (g). 


This completes the discussion on nontraditional features of the a-y scheme. In the 
following, it will be shown that this scheme can also be constructed from a completely 
different perspective. As a part of this construction, SE’s and CE’s of different types will 
be used and discussed. 
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In the new construction, the locations of mesh points (dots in Fig. 3(a)) are iden- 
tical to those used in the original construction. However, SE(j, n) is defined to be the 
interior of a rhombus centered at (j,n) (see Fig. 3(b)). CE(j, n) is the union of SE (j,n) 
and its boundary. Readers are warned not to confuse the sides of the rhombus with the 
characteristics of Eq. (2.22). Any one of these sides is simply a line segment joining two 
points of intersection (not marked by dots) of horizontal and vertical mesh lines. For any 
(x,t) 6 SE0», u(x,t) and h(x,t ), respectively, again are approximated by u*(r,t; j,n) 
and h*(x, t ;j, n ) which are defined by Eqs. (2.3) and (2.7), respectively. However, Eq. (2.5) 
will be derived from a consideration of flux conservation. 

Let Eq. (2.2) be approximated by 

/ h*-ds= 0 

Js(V) 

where V* is the union of any combination of CE’s. Because an SE is the interior of a CE, 
h * is not defined on S(V*), the boundary of V* . As a result, the above surface integration 
is to be carried out over a surface that is in the interior of V* and immediately adjacent 
to 5(V*). A necessary condition of Eq. (2.27) is that, for all (j, n) € 

® h* • ds — 0, 

J S(CE(j,n)) 

i.e., the total flux leaving any conservation element is zero. 

Note that the center of a current SE no longer sits on an interface separating two 
CE’s. It coincides with the center of a CE. Thus h* at one side of an interface is evaluated 
using information from one SE, while that at the other side is evaluated using information 
from another SE. As an example, h* at BC and B'C' depicted in Fig. 3(d), respectively, 
are evaluated using information from SE(j, n) and SE (j — 1/2, n — 1/2). Another necessary 
condition for Eq. (2.27) is the equality between the fluxes entering and leaving any interface. 
This can be seen by applying Eq. (2.27) separately to two neighboring CE’s, and then to 
their union. Obviously the local flux conservation relations at all interfacs, and within all 
CE’s (i.e., Eq. (2.28)) are equivalent to the global conservation relation Eq. (2.27). The 
equations representing the above conservation conditions are the numerical equations to 
be solved. Note that, in the current construction, a flux is not preassigned at an interface 
using an interpolation or extrapolation of information from both sides of this interface. 
The current method of interface flux evaluation obviously is different from that used in 
the finite volume method which was disscussed in Sec. 1. 

By using Eqs. (2.3) and (2.7), one concludes that, for any (x,t) 6 SE(/,n), the diver- 
gence of h* in £ 2 is 


(2.28) 


(2.27) 


def d[au*(x,t;j,n) — fj.du*(x,t]j,n)/dx] du*(x,t ;j,n) 

~ di + dt 


— a ( u x )" + («,)7 
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Because (u x )] and (u t )" are constants within an SE, Eq. (2.29) implies that V • h* is also 
a constant. Thus Eq. (2.28) coupled with Gauss’ divergence theorem implies that, within 

any SE, 

V • h* = 0. (2-30) 


Eq. (2.5) is a direct result of Eqs. (2.29) and (2.30). 

Note that Eq. (2.30) follows from Eq. (2.28) because u*(x, t ; j, n) defined in Eq. (2.3) is 
a first-order Taylor’s expansion. For a higher-order expansion, the condition that Eq. (2.30) 
being valid uniformly within an SE is stronger than Eq. (2.28). For the general case, the 
stronger condition should be imposed. Because Eq. (2.30) is the numerical analogue of 
Eq. (2.1), the imposition of the stronger condition ensures that, within an SE, the numerical 
solution uniformly satisfies the differential form of the conservation law Eq. (2.2). 

With the aid of Gauss’ divergence theorem, Eq. (2.30) implies that the surface inte- 
gration of h* over any closed surface located within any SE vanishes. As a result, 

(/ h * • ds = 0, and h* • ds = 0, (2.31) 

JS(AABC) JS(AA'B'C') 

where the triangles A ABC and A A'B'C' are those depicted in Fig. 3(d). Because the 
net flux of h* entering an interface from both sides vanishes, the sum of the flux leaving 
CE(j,n) through BC and that leaving CE(j - 1/2, n- 1/2) through B'C' vanishes. Thus, 
Eq. (2.31) implies that F_(j,n) = 0 where F-(j,n ) is defined in Eq. (2.11). Similarly, it 
can be shown that F+(j, n) = 0. 

Assuming Eqs. (2.3) and (2.7), it has been shown that both Eqs. (2.5) and (2.8) 
can be derived using Eq. (2.27). Conversely, Eq. (2.27) also follows from Eqs. (2.5) and 

(2.8) . Obviously both the forward marching scheme Eq. (2.14) and the backward marching 
scheme Eq. (2.22) can also be obtained by assuming Eqs. (2.3), (2.7), and (2.27). 

Note that the equivalence between Eq. (2.27) and the pair of equations Eqs. (2.5) and 

(2.8) hinges on the fact that V • h* = 0 within an SE of either type I or type II. As will be 
shown immediately, this condition can be used to simplify evaluation of the flux across a 
simple curve that lies entirely within an SE of either type. 

According to the top expression given in Eq. (2.29), V • h* =0 implies that there 
exists a function v(x.t :j,n) such that 


dv{x,t]j,n ) 

dt 


= au*(x,t ; j,n) — p 


du*(x,t;j,n) 

dx 


(2.32) 


and 


g^(r,f ;j,n) 

dx 


— ti (r , t , j , ti ) 


(2.33) 


for any (x,t) 6 SE (j,n). Substituting Eq. (2.6) into Eqs. (2.32) and (2.33), one concludes 
that, up to an arbitrary constant, 

i/>(x,t;j,n) = - - ~ g >) ~ a (t ~*”)] +2/i(f-t n )j (2.34) 

-uj[(x-xj)-a(t-t n )]. 
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Moreover, with the aid of Eq. (2.9), Eqs. (2.32) and (2.33) imply that 

g* ■ dr = dip. (2.35) 

Let (x,<) 6 SE(, 7 , n) and £ SE(j, n). Let T be a simple curve joining (x,t) and 

(x',t'), and lying entirely within SE(j, n) (see Fig. 4). Then Eqs. (2.10) and (2.35) imply 
that 

j h* ■ ds = ip(x',t'-, j,n) - ip(x,t (2.36) 

Here we assume that ds points to the right of T if one moves forward from (x, t) to (x',t') 
(see Fig. 4). Eq. (2.36) states that the flux of h* across the curve T is given by the difference 
in the values of tp at its two end-points. For this reason, ip(x,t;j,n ) will be referred to 
as the potential function associated with SE(j, n). Obviously, Eq. (2.12) can be obtained 
using Eq. (2.36). 

In [1], the a-fj, scheme is subjected to a thorough theoretical and numerical analy- 
sis on stability, dissipation, dispersion, consistency, truncation error, and accuracy. It 
is shown that it has many advantages over the MacCormack and the Leapfrog/ DuFort- 
Frankel schemes. Particularly, by using a discrete Fourier analysis, it is shown that the 
a-fi scheme is more accurate than the Leapfrog/DuFort-Frankel scheme by one order in 
both initial- value specification and the marching scheme itself. 

In conclusion, a model scheme has been constructed from two different perspectives 
using SE’s and CE’s of different types. Using either perspective, one can say that a numer- 
ical solution generated using the current framework satisfies (i) the differential form of the 
conservation law uniformly within an SE, and (ii) the integral form over any region that 
is the union of any combination of CE’s. In the author’s opinion, the second perspective 
that uses the SE’s and CE’s of type II depicted in Fig. 3 is more fundamental and thus was 
used in the initial development of the current method [1]. However, we believe that the 
first perspective is easier to use in constructing explicit schemes. Because most schemes 
to be constructed in the current paper are explicit, the first perspective will be adopted 
unless specified otherwise. 
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3. The a-e Scheme 


The inviscid a-ji scheme is neutrally stable and reversible in time. It is well known that 
a neutrally stable numerical analogue of Eq. (2.22) generally becomes unstable when it is 
extended to model the Euler equations. It is also obvious that a scheme that is reversible 
in time cann ot model a physical problem that is irreversible in time, e.g., an inviscid flow 
problem involving shocks. In this section, we assume p = 0 and attempt to modify the 
inviscid a-p scheme such that it can be extended to model the Euler equations. 

The current path of development is almost identical to that given in Sec. 2. We 
continue to assume Eqs. (2.3)— (2.7), and use SE’s of type I depicted in Fig. 2. In addition 
to // = 0, the only other modification is the replacement of the assumption F±(j, n) = 0 


by 


F±(j,n ) = ± — 


i/ 2 )(ax ) 2 

_ 


(<*«,)?, 


(3.1) 


where e is a parameter independent of numerical variables, and 

<*.)? = ( 1 / 2 ) [<-)&£ + («.)?:$] - («;;# - «;:#) /**• (3-2) 


In other words, we add two terms of the same magnitude but with opposite signs, re- 
spectively, to the right sides of the original conservation conditions F+(j,n ) = 0 and 
F-(j,n ) = 0. The beauty of this modification will be fully explained later in this section. 
For now it s uffi ces to say that this modification injects a higher-order finite-difference error 
into the inviscid a-p scheme. It breaks the space-time symmetry of the latter. In turn, 
numerical diffusion is introduced as a result of this symmetry breaking. Because the mag- 
nitude of the terms added in this modification is controlled by e, numerical diffusion is 
controlled by e in the modified scheme just as physical diffusion is controlled by p in the 
a- (i scheme. Note that, as a result of Eq. (3.1) and the assumption p = 0, the modified 
scheme is characterized by two parameters a and e. Thus, hereafter it will be referred to as 
the a-e scheme. Also note that, because there is no upwind bias in the a-e scheme, upwind 
bias is not the source of numerical diffusion. Additional remarks on Eqs. (3.1) and (3.2) 


are: 

(a) By definition, F+(j,n) and FL(j, n) represent total fluxes leaving CE+(j, n) and 
CE -(j,n), respectively (see Figs. 2(c) and 2(d)). Because F±(j, n) 7^ 0 if e 7^ 0, 
CE + (j, n) and CE_(j, n) generally axe no longer conservation elements in the a-e 
scheme. 

(b) Let CE(j,n) be the union of CE+(j,n) and CE_(j, n) (see Fig. 5(b)). Note that this 
definition of CE(j,n ) differs from that given in Sec. 2 and depicted in Fig. 3(c). Let 

F(j, n) d = <f h* ■ ds. (3-3) 

Js{CE(j,n )) 

Because the net flux entering the interface separating CE+(j, n) and CE_(j, n) is zero. 
F(j,n) is the sum of F+(j, n) and F_(ji,n). With the aid of Eq. (3.1), we have 

F(j,n) = F+(j,n) + F-(j, n) = 0, (3.4) 


15 


i.e., the total flux leaving CE(/, n) vanishes. As a result, CE(j, n) is a conservation 
element in the a-e scheme . Note that Eq. (3.4) leads to a global conservation relation 
in the form of Eq. (2.27) where V* is the union of any combination of these new CE’s. 
(c) Because £ = 0 if p = 0, Eq. (3.4) coupled with Eq. (2.12) implies that 


«” = (1/2) [(1 + ■'K-l'/’ + 0 - "K 


n — 1/2 
>+ 1/2 


Ax(l — l/ 2 ) 

8 


[(“*)_,•_ 1/2 ( u x )j 


n-l/2\ 


U *Jj+ 1/2 


(3.5) 


Thus, u” is independent of e. 

(d) Because (u x )J±i/ 2 i s a numerical analogue of du/dx at (j ± 1/2, n — 1/2), the simple 
average 


( 1 / 2 ) [M";l ' 2 



is a numerical analogue of du/dx at (j, n — 1/2), the midpoint of a line segment joining 
(j + 1/2, n - 1/2) and (j - 1/2, n - 1/2) (see Fig. 2(a)). Note that (j,n - 1/2) ^ ft if 
(/, n) e ft. Also note that 


( U 


n— 1/2 
. 7 + 1/2 


— U 


n— 1/2 
J-l/2 



AX 


is a central-difference analogue of du/dx at (j,n — 1/2). Thus, (du x )j represents the 
difference of two numerical analogues of du/dx at the same mesh point (j, n — 1/2). 
By using Taylor’s expansion at (j, n — 1/2), it can be shown that (du r )" = O [(ax) 2 ] , if 

( u x)j±i /2 are identified with du(xj± 1 / 2 ,t n ~ 1 ^‘ 2 )/dx, respectively. Hereafter a quantity 

is denoted by O [( ax)^| if there exists a constant C > 0 such that the absolute value of 
this quantity < C |ax| ? for all sufficiently small |ax|. Note that we have constructed an 
expression of 0[(ax) 2 ] without explicitly introducing the factor (ax) 2 . This natural 
construction leads to the simple stability conditions to be given in Eq. (3.14). It is 
possible only because there axe two discrete variables u” and (u x )J associated with 
the mesh point (j.n). 

(e) Eq. (3.1) could have been written as F±(j, n) = ±e'(du x )” with e' = e(l — i/ 2 )(ax) 2 /4. 
However, this simplified expression would lead to much more complicated equations 
later. 


This completes the discussion of Eqs. (3.1) and (3.2). Now, let 1 — v 2 ^ 0. Then 
Eqs. (2.12), (3.1), and (3.2) can be used to obtain the current counterparts of Eqs. (2.14) 
and (2.18). They are 


q(j,n) = A/+ q(j - 1/2, n - 1/2) + M_ q(j + 1/2, n - 1/2) (1 - i^ 2 ± 0) (3.6) 


and 


JO, n + 1) = (M+) 2 g(j - 1, n) + (M+M- + M_.M+)q(j,n)+(M-) 2 q(j + 1, n) 

( 1 -^ 0 ), 


(3.7) 
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respectively. Here 


M + *= ( 1 / 2 ) 


and 


M- = (1/2) 


1 + 1/ 

l-i/ 2 ^ 

e — 1 

2e — 1 + v ) 

l-i/ 

-(1 - v 2 ) \ 

1-6 

2e - 1 - v ) 


(3.8) 


(3.9) 


Obviously, M± = Q± if e = 0 and ( = 0. Furthermore, the limiting condition given 
in Eq. (2.19) is still valid if we assume that e = e(At) and lim^t— o e(Af) = 0. However, 
unlike the a-p scheme, the a-e scheme is not a two-way marching scheme if e ± 0. 

Eq. (3.6) represents a pair of equations. The first is Eq. (3.5). With the aid of 
Eqs. (2.5) and (2.13), the second equation can be expressed as 

(“•)? = (“£■/* - “i-1/0 /ai + (26- !)(*.,)? . (3.10) 


«&./, = « *«£ + (A< /2)(u, r^ljl (3.11) 

i.e., a first-order Taylor’s approximation of u at (j ± 1/2, n). Thus, the expression 

on the right side of Eq. (3.10) is the sum of a central-difference approximation of du/dx 
at (/, n) and the extra term (2e - l)(du x )y. Because (d«x)” = 0[(ax) 2 ], the presence of 
this extra term will not lower the order of accuracy of the entire sum as an approximation 
of du/dx at ( j,n ). Also note that this extra term vanishes when e = 1/2 while the term 
associated with (du x )” in Eq. (3.1) vanishes when e = 0. 

Next we shall study the influence of e on the stability and numerical diffusion of the 
a-e scheme. Let G^} and be the principal and spurious amplification factors of the 
a-e scheme, respectively. Then, it will be shown in Sec. 6 that 



[A±(e,t',«)l 2 , 


(3.12) 


with 

A±(e, i/, 0) *== e cos(^/2) — iu sin(0/2) + a/( 1 — e) [(1 — e)cos 2 (#/2) + (1 — i/ 2 )sin (0/2)] • 

(3-13) 

Here 0 , —x < 0 < 7T [1, p.30], is the phase angle variation per Aar. Also it will be proved 
that 

0 < e < 1, and v 2 < 1 (314) 
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are necessary conditions for the stability of the a-e scheme. Thus, Eq. (3.14) will be 
assumed in the remainder of this section. 

It was pointed out in Sec. 2 that the amplification factors of the Leapfrog scheme are 

identical to those of the inviscid a-/z scheme. Because the latter scheme is a special case of 

( 2 ) . ^ 
the a-e scheme with e = 0, become the amplification factors of the Leapfrog scheme 

when e = 0. This fact can be reverified by comparing Eqs. (2.21), (3.12) and (3.13) with 
£ — 0 and e = 0. 

Also, we have 

= cos(0/2) — ivsin(0/2). (3.15) 

Thus, G\ = G_ when e = 1 . Moreover, it is shown in Appendix A that the coalesced 
amplification factor is identical to that of the Lax scheme. Note that, like the Leapfrog 
scheme, a solution of the Lax scheme is also composed of two decoupled solutions with 
each being associated with a mesh that is staggered in time. However, because the Lax 
scheme is a two-level scheme, it does not have a spurious amplification factor. 

Thus, at one extreme, i.e., when e = 0, G±* become the amplification factors of the 
Leapfrog scheme, which is free of numerical diffusion. At another extreme, i.e., when e = 1, 
GY and G_ ’ coalesce into one and it becomes the amplification factor of the Lax scheme, 
which is notorious for its large diffusive errors. Prom the above observations, one may infer 
the following conclusion that will be established shortly, i.e., the a-e scheme becomes more 
diffusive as the value of e increases. Note that, because the Lax scheme is very diffusive 
and uses a mesh that is staggered in time, a two-level scheme using such a mesh is usually 
associated with a highly diffusive scheme [18]. The a-e scheme demonstrates that it can 
also be a scheme with no diffusive error! 

As a result of Eq. (3.14), the expression under the radical sign in Eq. (3.13) is non- 
negative. Thus, it can be shown that 


1 - |G± J | = X±(€,i/,0) *= ej(l - i/ 2 )sin 2 (0/2) + 2cos(0/2)x 

(3-16) 

(1 - e) cos(0/2) =F yj{l - 7 ) [(1 - e)cos 2 (0/2) + (1 - i/ 2 )sin 2 (0/2)]| 1 


Because solutions to the physical equation Eq. (2.22) do not dissipate with time, a numer- 
ical analogue to Eq. (2.22) is said to be free of numerical diffusion if its solutions also do 
not dissipate with time, i.e., its amplification factors are of unit magnitude. As a result, 
numerical diffusion of the a-e scheme may be measured by 1 — |G±*|, i.e., x±(e, v, 0)- Ob- 
viously the a-e scheme is free of numerical diffusion if e = 0. Also, by using Eqs. (3.14) 
and (3.16), it is shown in Sec. 6 that, for all 6 with —tt < 9 < it, and all e and v satisfying 
Eq. (3.14), we have 

0 < X+(e, v,9) + 4e(l - e)cos 2 (^/2) < x_(e,i/,0) < min{l,4e}, (3.17) 
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and 


0 < x+(e, v, 6) < e (1 - v 2 ) sin 2 (0/2). 


(3.18) 


The significance of Eqs. (3.17) and (3.18) is discussed in the following remarks: 

(a) Because 0 < x±(t,v, e ), Eq. (3-16) implies that |G ( ± } | < 1. It is proved in Sec. 6 
that this result and other considerations lead to the conclusion that Eq. (3.14) is also 
sufficient for stability. 

(b) For a numerical analogue of Eq. (2.22) that has both principal and spurious ampli- 

fication factors, a numerical solution with periodic boundary conditions is the sum 
of a principal solution and a spurious solution [1, p.37]. Only the principal solution 
contributes to the accuracy of the scheme. Given a smooth initial condition, the spu- 
rious solution at t = 0 generally is very small compared with the principal solution. 
Also, the behaviors of the principal and the spurious solutions as functions of time are 
determined by the principal and spurious amplification factors, respectively. Because 
0 < e(l — e) if 0 < e < 1, Eq. (3.17) implies that x+( e > v ^) ^ Thus, the 

spurious solution will dissipate not slower than the principal solution. Let e be not 
too close to 0 or 1. Then Eq. (3.17) also implies that the Fourier components of the 
spurious solution with smaller |0|, i.e., longer wavelength, will dissipate much faster 
than those of the principal solution. In other words, the spurious solution will rapidly 
disappear from the long-wavelength components of a numerical solution. Note that 
X_(l/2, i/,0) = 1. Thus, the long-wavelength components of the spurious solution are 
annihil ated almost completely in a single time step if e = 1/2, i.e., if the last term in 
Eq. (3.10) is dropped. 

(c) The upper bound of x+(«»*'»0) given in Eq. (3.18) is proportional to sin 2 (0/2). As 
a result, the long-wavelength Fourier components in the principal solution are nearly 
free of numerical diffusion. On the other hand, short-wavelength components may 
decay rapidly. 

(d) For a fixed e, Eq. (3.18) implies that the principal solution is more diffusive for a 
smaller |i/|. To compensate, one may choose a small e for a small \v\. One may even 
choose e to be a monotonic function of |i/|, subjected to the condition 0 < € < 1. 

(e) Eqs. (3.17) and (3.18) imply that, for all v with v 2 < 1 and all 6 with -n < 0 < tt, 
we have 


0 < x+(e,M) < and 0 < x-(«, v, 0) < min{l, 4e}, (3.19) 

which, according to Eq. (3.16), is equivalent to 

1 — e<|G^|<l, and 1 - min{l,4e} < ) | < 1. (3.20) 

As a result, by choosing t small enough, both |G^| and |G^| can be confined within 
an arbitrarily narrow range. As noted previously, the spurious part of a numerical 
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solution generally is insignificantly small assuming a smooth initial condition. It 
does not contribute to accuracy and usually dissipates faster than the principal part. 
Thus, our primary concerns is how the principal part dissipates. From Eq.(3.20), one 

concludes that, for any e with 0 < e < 1, |G+ | will be bounded uniformly from below 
by a positive number 1 - e for all 1 / with v 2 <1 and all 8 with -x < 8 < tt. By choosing 
an e of proper magnitude, one can suppress artificial numerical oscillations without 
causing large diffusive errors for any combination of v and 0. This fact contrasts 
sharply with what one expects from typical classical schemes which are usually very 
diffusive with respect to certain v and 0, while not at all with respect to other v and 
0. As an example, we consider the Lax-WendrofF scheme [4, p. 101] . Its amplification 
factor is of unit magnitude, for all 8 at v = 0, or v — 1. On the other hand, the 
amplification factor = 0 if v 2 — 1/2 and 0 = ir. 

In nonlinear flow solutions, e.g., shock-tube solutions to be discussed in Sec. 8, 
analogues of v axe dependent on local velocity components. Thus, they may vary 
from one location to another. Also, at some neighborhood, the Fourier spectrum of 
the local solution may have peaks spread over a wide range of 8. Thus, for a numerical 
analogue of Eq. (2.22), a large variation in numerical diffusivity with respect to 0 and 
v generally means that numerical solutions obtained using its nonlinear extensions 
will suffer annihilations of sharply different degrees at different locations and different 
8. Such selective annihilations may cause large distortions of numerical solutions [19], 

This completes the discussion of Eqs. (3.17) and (3.18). In conclusion, the a-e scheme 
has been constructed to solve Eq. (2.22). It has the unique property that numerical diffu- 
sion can be controlled by a parameter e. Because neither characteristics- based tec hni ques 
nor knowledge about the upwind direction is used in the construction of the a-e scheme, as 
will be shown in the next section, it can be easily extended to model the Euler equations. 
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4. The Euler Solver 


We consider a dimensionless form of the 1-D unsteady Euler equations of a perfect 
gas. Let p, v, p, and 7 be the mass density, velocity, static pressure, and constant specific 
heat ratio, respectively. Let 


m = p, u 2 = pv, u 3 =p/(7-l) + (l/2)pt> 2 , 

fi = u 2 , 

h = (7 ~ 1)«3 + (l/2)(3 - 7 )(u 2 ) 2 /ui, 


and 


/ 3 = 7U2U3/W1 - (l/2)(7 “ 1)(«2) 3 /(«i) 2 - 
Then the Euler equations can be expressed as 

du m df n 


+ -4r- = 0 , 


dt dx 

The integral form of Eq. (4.5) in space-time E 2 is 


/ 

JS(V) 


h m ' ds — 0, 


m = 1 , 2 , 3 . 


m = 1,2,3, 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 


where h m = (/ m , u m ), m = 1,2,3, are the space-time mass, momentum, and energy 
current density vectors, respectively. 

As a preliminary, let 

fm,k d = dfm/du k , m,k = 1,2,3. (4.7) 

Let F be the matrix formed by m,k = 1,2,3. Then Eqs. (4.2)-(4.4) imply that 

/ 0 


F = 


(3-7) f U2 

. U 1 


(7-1) ( Hi ) _7^H1 

17 7 (u,) 2 


1 


0 

(3-7)- 

u 1 


7-1 

u 3 3(7 - 1 ) fu 

7 Ul 2 \v 


u 2 
7— 


(4.8) 


Let c be the sonic speed. Then 


c = \ 

7(7 - 1) 

u 3 1 / 1 

\ 2 " 
-) 



ui 2 \ti 

‘i ) 


(4.9) 
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Let G be the 3x3 matrix defined by 


/ 


G d = 


u 2 
« 1 


KS ! 


(« 2) 2 


J*l_ 

s/2c 

u 2 _ Ul_ 

V^c 


U2 -j_ u l c 


(U 2 ) 2 


«1 

\f2c 


«2 «1_ 
V^c + vS 


U 2 

+ -4 + 


\ 


ItjC 


2\^CU! V5 v/2(7 - 1) 2V2CU! y/2 y/2(y - 1) / 

Then the inverse of G is given by 


(4.10) 


G = 


1 (7 “ 1 

(aV 

(7 - l)u 2 

(7 — 1) \ 

2 c 2 1 

Ui/ 

C 2 Ui 

c 2 

(7-1)(u 2 ) 2 

u 2 

1 (7 - l)u 2 

7-1 

2 ^c ( Ul ) 3 

V^i ) 2 

\/ 2 ui \/ 2 c(tti ) 2 

s/2c Ul 

(7 — l )( u 2 ) 2 

u 2 

1 (7 - l)u 2 

7-1 

V 2\/2c(ui) 3 

V2( Ul y 

V^ui V^c(u !) 2 



(4.11) 


For any numbers ai, 02 , . . o n , let diag(ai,a 2 , . . . , a n ) denote the diagonal matrix with , 
02 , . . ., a n being the diagonal elements on the first, second, . . ., and n-th rows, respectively. 
Then, by using Eqs. (4.8)-(4.11) and v = u 2 /ui, one has 


G 1 F G = diag(v , v — c, v 4 - c). 


(4.12) 


Consider SE’s of type I depicted in Fig. 2 . For any ( x,t ) € SE(j,n), u m (x,t), 

and h m (x,t ) are approximated by u^(x,<;j,n), /^(x,t ;j,n), and h* m (x,t-, j,n), respec- 
tively. They will be defined shortly. Let 

def 

= {u m )j +(u mx )”(x - Xj) + (u m t)j(t -t n ), m = 1,2,3, (4.13) 

where (u m )", (u mi )", and (u mt )" are constants in SE(j,n). Obviously, they can be con- 
sidered as the numerical analogues of the values of u m , du m /dx , and du m jdt at (xj,t n ), 
respectively. 

Let (/ m )” and (/ m ,jt)> denote the values of f m and f m ,k, respectively, when u m , 
m = 1,2,3, respectively, assume the values of (u m )", m = 1,2,3. Let 

3 

(fmx)j — ^ ^ {fm,k)j (ukx)j , m = 1,2,3, (4.14) 

fc=l 
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(/-.)? =' £(/■»,*)”(“ «)?, 


m = 1,2,3. 


(4.15) 


d/m \T^ r 

2 -^ •' m,< 


(4.16) 


d/m _ r dufc 

dt - 2-j Jm ' k dt ’ 


(4.17) 


(fmz)j (fmt)j can be considered as the numerical analogues of the values of df m /dx 
and df m /dt at (xj,t n ), respectively. As a result, we assume that 

/* (x,t;i,n) = (/„)" + (/ mI )”(x - Xj) + (/„,)"(( - ("), m = 1,2,3. (4.18) 


Because h m = (/ m , u m ), we also assume that 

hm^Xft j j , w) = (/ m (x,t,J, Tl), U m (x,t , *i)), 


m = 1,2, 3. (419) 


Note that, by their definitions, (i) ( f m )] and ( , m = 1,2,3, are functions of (u m )?, 
m = 1,2,3, (ii) ( f mx )”, m = 1,2,3, are functions of (u m )” and (u mi )”, m = 1,2,3, and 
(iii) (/mt)” are functions of (u m )" and (u m <)y , m = 1,2,3. 

Moreover we assume that, for any (x,t) £ SE(j,n), u m = u^(x,t ;j,n) and / m = 
fm(x,t;j,n) satisfy Eq. (4.5), i.e., 


du^(x,t;j,n) d/^(x,f ;j,n) _ 
dt dx 


According to Eqs. (4.13) and (4.18), Eq. (4.20) is equivalent to 

(‘ «ml )? = “(/-)”• 


(4.20) 


(4.21) 


Because (f mx )] are functions of (u m )" and (u mr )", Eq. (4.21) implies that (u mt )] are also 
functions of (u m )” and (ti mi )". From this result and the facts stated following Eq. (4.19), 
one concludes that the only independent discrete variables needed to be solved in the 
current marching scheme are (u m )j and (u mi )”. 

From Eq. (4.20), one concludes that the generalization of the potential function 
ifi(x,t ; j,n) introduced in Sec. 2 to the current solver are xl> m (x,t;j,n ), m = 1,2,3, which 
satisfy 

= (4.22) 


23 



fa = u m OM;j,n). 


(4.23) 


Substituting Eqs. (4.13) and (4.18) into Eqs. (4.22) and (4.23), and using Eq. (4.21), one 
concludes that, up to an arbitrary constant, 


*m(x,t-,j,n) = (/„)"(< - n - («»)?(* - *>) + (1/2 )(/„«)“(( - <") 2 

- (l/2)(«„,)J(x - x;) 2 + (/«,)?(* - x;)(f - <”)• 

By using an argument similar to that leading to Eq. (2.36), one concludes that 


(4.24) 


j h" m -ds= 4’ m (x\t';j,n) - r/> m (x,t ;j,n). 


(4.25) 


Here T is a simple curve joining (x,t) and (x',f'), and lying entirely within SE(j,n). We 
also assume that ds points to the right of T if one moves forward from (x,f) to (x',f'). 

As in the a-e scheme, we assume that the flux of is conserved over CE(j, n), i.e., 


' S(CE(j,n)) 


Combining Eqs. (4.25) and (4.26), one has 


h* m • ds = 0. 


(4.26) 


*l>m(xj - ax/ 2 ,t n ; j,n ) - ip m (xj + ax/ 2, t”; j,n ) 

+ ipm(xj- 1/2 + ax/ 2 ,< n_1/2 ; j - 1/2, n — 1/2) 

- tM*j-i/2,< n_1/2 + At/2; j - 1/2, n - 1/2) (4.27) 

+ Vi m (x j+ i /2 ,t n ~ 1/2 + At/2; j + 1/2, n - 1/2) 

- 0 -m0r;+i/2 - a x/2,t n ~ 1/2 ; j + 1/2, n - 1/2) = 0. 

Substitution of Eq. (4.24) into Eq. (4.27) yields 

( u m)j = ^ ( Um )>-l/2 + ( Sm )j-\/2 (-Sm ) ] > (4.28) 

where, for all (j, n) G Q, 

(*-)J = + + ”> = 1,2,3. (4.29) 

Eq. (4.28) forms the first half of the current marching scheme. The second half which 
solves (u mz )] will come from a generalization of Eq. (3.10). 
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For all (j, n) G ft, let 


/ j„ \n def ^ 

(aU mx )j — 


2 L 


(u 


,n-l/2 / 

mx)j+\/2 * \ u t 


,)£$] - [(x™)”;, 1 / 2 - M’ll'/l] A*. (4.30) 


U4.1VJ. y 

K.W = («-)£# + ( i( /2)(«mi)"±i/2 1 (431) 

for m = 1,2,3. Because Eqs. (4.30) and (4.31) are the generalizations of Eqs. (3.2) and 
(3.11), respectively, a natural generalization of Eq. (3.10) is 

(«-«)? = [(0*1/2 - /« + (2' - !)(<**-.)“. m = 1. 2. 3, (4.32) 

where c is a parameter independent of numerical variables. Note that the last term 
in Eq. (4.32) vanishes if e = 1/2. The marching scheme presented in [2] is formed by 
Eqs. (4.28) and (4.32) with e = 1/2. 

To construct a larger class of generalizations to Eq. (3.10), for all (j, n) G ft, let 


(*-)* = \ 


^UmTj+X/l + (“m)"_ 1 1 /2 J 


- 1 / 2 ] 


m = 1,2, 3. 


(4.33) 


Let (e m )j, m = 1,2,3, be parameters that can be functions of (u m )", m — 1,2,3. There 
can be many choices of these functions. Let ( Qmk )j be the value of the (m, k )— element of 
the matrix G when tt m , m = 1,2,3, respectively, assume the values of m = 1,2,3. 

Similarly, let (<7~]t)” be the value of the (m, fc)-element of the matrix G -1 when u m , 
m = 1,2,3, respectively, assume the values of (u m )”, m = 1,2,3. Let 


(e-OJ = m ’ k = 1>2 ’ 3 ' 


(4.34) 


i=l 


Then Eq. (3.10) can be generalized as 


(Umx)j = f( u m)>+l/2 ~ ( U m)>- 1/2 / A;r + [2(^mi)j ( du kx )j , (4.35) 

L fc=l 


where m = 1,2,3, and <5 m jt is the kronecker-delta symbol. 

Consider the special case in which, for all (j, n) G ft, (^i)j == (<h)j = (^)j • Let 
(« m )” = (e)", m = 1,2,3. Then (e mJt )” = {e)]6 mk , and thus Eq. (4.35) is reduced to 


(Umx)” = ( u m)j+l/2 - ( U m)> — 1/2 /ax + [2(c)" — l] (du mz )", TU — 1,2,3, 


(4.36) 


Note that Eq. (4.36) reduces to Eq. (4.32) if (e)J = e for all (j,n) G ft. 
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The current Euler solver is a straightforward extension of the a-e scheme. Eqs. (4.28) 
and (4.35), which form the marching scheme, will be converted into a matrix form similar 
to Eq. (3.6). Eqs. (4.43) and (4.44), to be obtained during this conversion, will also be 
used in a stabilty analysis presented in Sec. 6. To proceed, note that f m , m = 1,2,3, are 
homogeneous functions of degree 1 [20, p.ll] in the variables u m , m = 1,2,3. Thus 


(/-)* = £(/■»,*)?(«*)?• 

k=l 

Also, by using Eqs. (4.15), (4.21) and (4.14), one has 

(/-.)? = - E »)?(«■*.)?• 

k= i ;=i 

Substituting Eqs. (4.37) and (4.38) into Eq. (4.29), and using the definitions 

«.)”=' m = 1,2,3, 

and 

(O* = £;(/».*)■. m,k— 1,2,3, 

one has 


(4.37) 


(4.38) 

(4.39) 

(4.40) 


(*-)" = E(£, *)’(«*)? + E 


*=1 


Jk=l 


- E(£,<>"(tf»)" 


/=1 


(<&)". m = 1,2,3. (4.41) 


Substituting Eqs. (4.21) and (4.14) into Eq. (4.31) and using Eqs. (4.39) and (4.40), one 
also has 

(•4W = («.)£$ - 2 *■ < 4 - 42 ) 

t=i 

With the aid of Eqs. (4.39)-(4.42) and (4.30), Eqs. (4.28) and (4.35) can be expressed 
as 


(«m)” — g ^ mk + 1/2 ( U *)j- 1/2 


*=1 


+ 




i=l 


,+ '\ n_ 1/2 
- 1/2 


/ + ] 

( U kx )>-1 


+ [*m* - (/m.^i+Z/z 2 ] («*)". 


i-l/2 

^+ 1/2 


i=i 


( u + -l”-!/ 2 ! 

\ U kz)j+l/2 J 5 


(4.43) 


m = 1 , 2 , 3 , 
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and 


= 5 e{ [(«-*)■ - M [(*■*);:$ 

*=1 l 


(«fc) 


n— 1/2 

i+1/2 


+ 

+ 


- 1/2 
/2 


[2(e ml )- - *». + (/£»),":$] wf,);:,' 

[2(e mt )" - S mi - j- 


m = 1,2, 3, 


(4.44) 


respectively. 

For all (j, n ) G ft, let u” and (u+ )", respectively, denote the column matrices formed 
by (« m )" and m = 1,2,3. Let (£)” and (F + )”, respectively, denote the 3x3 

matrices formed by (e m k)j and (/* k )j, m,k = 1,2,3. Let I denote the 3x3 identity 
matrix, i.e., the matrix formed by m,k — 1,2,3. For all (j,n) G ft, let 


q (j, n) 


def 



(4.45) 


M + (j,r,) d = (1/2) 


/ / + (B + )”:$ 


and 


M- (/, n) d = ( 1/2) 


I - (E)] 


- \f+)X, 


- 1/2 
/2 


T 2 


\ 


V (£)? - / 2(B)? -1 + [F+T-'J/l ) 


(4.46) 


[(f + )” 


n— 1/2 

i+1/2 


i 2 


/ \ 


2(E)? - / - (F+)];;/l / 


(4.47) 


By their definitions, q(j,n) is a 6 x 1 column matrix while M+ (j, n) and M_ (j, n) are 
6x6 matrices. With the aid of Eqs. (4.45)-(4.47), Eqs. (4.43) and (4.44) can be cast into 
the matrix form 


q (j, n) = M+ (j,n)q(; - 1/2, n - 1/2) + M_ (j,n)q(j + 1/2, n - 1/2). (4.48) 

Note that the set of equations given in Eqs. (4.45)-(4.48) and the set given in Eqs. (2.15), 
(3.8), (3.9) and (3.6) are very similar in their forms. 
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Recall that both v and c are functions of u m , m — 1,2,3. For all SE (j, n), let vj and 
c", respectively, denote the values of v and c when u m , m = 1,2,3, respectively, assume 
the values of ( u m )j , m = 1,2,3. It will be shown in Sec. 6 that the marching scheme 
defined by Eq. (4.48) can be linearized and decoupled into three pairs of equations with 
each pair being in the form of Eq. (3.6). This decoupling and other considerations lead to 
the conclusion that Eq. (4.48) is stable if, for all (j, n) £ ft, 

(*«.«)"<!, and 0<(e m )?<l, m = 1,2,3, (4.49) 


where 

(*«.,); =' (i*;i + ia?i) (4-50) 

We conclude this section by introducing some possible modifications to the above 
solver. Note that (u' m )" ±1 / 2 , by its definition, represents a finite-difference approximation 
of u m at ( j ± 1/2, n). As a result, 

«,)" = [(OJh/j - (<C)?-,/ 2 ] /**, ™ = 1.2.3, (4.51) 


respectively, are the central-difference approximations for du m /dx, m = 1,2,3, at (j, n). 
Note that (u£, x )" is the first term on the right side of each of Eqs. (4.32), (4.35) and (4.36). 
The above central-difference approximation is valid as long as no discontinuity of u m (or 
its derivatives) occurs between ( j - 1/2, n) and (j + 1/2, n) (see Fig. 5). In the following 
discussion, we develop alternates which are valid even in the presence of discontinuity. 


Let 


(«m.±)? = ± 


„ def ( U m)j±l/2 ( U ™)> 


ax/2 


m = 1, 2, 3, 


(4.52) 


where (« m )y can be obtained from Eq. (4.28). Because ( u m)7-i/2’ ( u m)yfi/ 2 ’ 

are the numerical analogues of u m at (j — 1/2, n), (j, n) cind (j -1- 1/2, n), respectively, 
(u mi _)" and (u mI+ )" are two numerical analogues of the value of du m /dx at (j,n) with 
one being evaluated from the left and another from the right. Note that 


( u mx)j ~ o [( u mx-)" + («mr+)"] • 


(4.53) 


In case a discontinuity occurs between (j, n) and (j + 1/2, n) but not between (j,n) 
and (j - 1/2, n), one would expect that |(u mt+ )”| > |(u mx _)’?|. Moreover, because (j,n) 
and (j — 1 /2, n) are on the same side of the discontinuity while (j, n ) and ( j + 1/2 , n) are on 
the opposite sides, (u m x)" should be a weighted average of (u mr+ )" and (u mx _)" biased 
toward the one with the smaller magnitude. 

As a result of the above considerations, (u^ r )" can be replaced by 

(t c*)" d = Wo ((ti m ,-)J,(ti mr+ )7;a) , m = 1,2,3. (4.54) 
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Here a is an 

(ii) 


adjustable constant and the function W 0 is defined by (i) W^ o (0, 0, a) — 0 and 


tF 0 (x_,x + ;a) 


|X-f |°X— + |g-| a x + 
|x + |° + |x_|“ 


(|x+| + |x_| > 0) 


(4.55) 


where x+ and x_ axe any two real variables. Note that W 0 (x-,x+;a) — (x_ + x+)/2, 
i.e., (uj*)* = (u c mx )], if a = 0 or |x_| = |x + |. Also the expression on the right 
side of Eq. (4.55) represents a weighted average of x_ and x+ with the weight factors 
|x + | a /(|x + |° f + |x_|") and |x_| a /(|x + | a + |x_|“). Fora > 0, this average is biased toward 
the one am ong x+ and x_ with the smaller magnitude. For the same value of |x+| and 
|x_|, the bias increases as a increases. Thus, we should always choose a > 0. 

Note that the special weighted averages Wo(x_,x+;l) and W 0 (x_,x+; 2) are used in 
the slope-limiters proposed by van Leer [21] and van Albada [22], respectively. 

The above modification, i.e., (u c mx )J replaced by («*•*)?, is first given in [2]. It is 
shown in [2] and also Sec. 8 of the current paper that it is an efficient tool to suppress 
overshoots and/or numerical oscillations near a discontinuity. Moreover, because ( u mx ±) r - 
axe constructed using only the data associated with the mesh points (j 1/2, n 1/2) and 
(j + 1/2, n - 1/2), the effect of this modification is highly local, i.e., it generally will not 
cause the smearing of shock discontinuities. 

However, there may be a price to pay for the above modification. Because a fractional 
power is costly to evaluate, so is W 0 (x_,x+;a) if a is not an integer. Moreover, because 
the bias of this weighted average increases with a , a situation may arise such that the use 
of an a with |a| < 1 may be desirable. To obtain a computationally efficient weighted 
average of arbitrary small bias, let 


W(x_, x+; a, (3) = (1 - 0)W o (x 


x + ;0) + 0W o (x-,x+; a), 


(4.56) 


where /? > 0 is an adjustable weight factor, and a generally is an integer. Because 
W / £) (x_,x+; 0) is the simple average of x_ and x + , Eq. (4.56) defines a linear weighted 
average of this simple average and the nonlinear weighted average defined in Eq. (4.55). 
Obviously, W(x_,x + ;a,/?) = (l/2)(x_ +x + ) if x_ = x + . Furthermore, because 


TF 0 (x_,x + ; -a) 


|x + [ a x+ + |x-| a x_ 

|*+|° + \x-\ a 


(|x+| + k-l > 0), 


(4.57) 


alternatively, W(x_, x + ; a, ft) can also be expressed as 

W(x-,x+;aJ ) = W 0 {x-,x+,-,a)+ W 0 (x_, x + ; -a). (4.58) 

The application of the more general modification, i.e., (u^i )> replaced by 


«x)” d = ^((“m*-)",(«m*+ )";«,/?), m = 1, 2, 3, (4.59) 

will be demonstrated in Sec. 8. 

Finally, note that W(x_, x + ; a, 0) can be further generalized by a linear weighted 
average of several W 0 (x_, x + ; a) with different values of a. 
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5. The Navier— Stokes Solver 


We consider a dimensionless form of the 1-D unsteady Navier-Stokes equations of 
a perfect gas [4, pp. 191-193]. (Note: the expressions on the right sides of the last three 
equations in Eq. (5-47) of [4] have incorrect signs in the earlier versions. The conduction 
heat-flux vector should be proportional to the negative of the gradient of temperature.) 
These equations are extensions of the Euler equations defined in Sec. 4. Thus, unless 
specified otherwise, the symbols, definitions, and equations given there will be used in this 
section. 

Let Rei and Pr denote the Reynolds number and Prandtl number, respectively. They 
axe assumed to be nonnegative constants. Let 


fi = 0, 


7 def 
J2 = 


« 2 


3 Rei u\ ' 


and 


7 de{ 
/ 3 = 


ZR&L 


(S : 


+ 


Rei Pr 




(«2> 


2 1 


Ui 2( U1 )2J 


Then, the Navier-Stokes euations can be expressed as 

d 2 f n 


du m dfn 


= 0 , 


m = 1,2,3. 


dt dx dx 2 
The integral form of Eq. (5.4) in space-time E 2 is Eq. (4.6) with 

m = 1, 2, 3. 


h m d = (f m - df m /dx, u m ), 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 


As a preliminary, let 


fm,k d = dfm/dllk, 


m,k = 1,2,3, 


(5.6) 


and 

def 4 def 7 , def /r 

Tl = 3 RT l ' T2= R^Tr' “ d TI = T2 - T - < 5 ' 7 ) 

Let F denote the 3x3 matrix formed by fm,ki m,k = 1,2,3. Then Eqs. (5.1)-(5.3) imply 
that 


( 


F = 


(“2) : 

v T3 (^ 


0 

0 

0 
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Again we consider SE’s of type I depicted in Fig. 2. For any (x, t) (E SE(j, n), 
u m (x,t), / m (x,f), / m (x,t), and h m (x,t ), respectively, are approximated by u* m (x,t;j,n), 
fn(x,t-,j, n), /^(x, t ; j, n), and C(x,<;j, n). u* m (x,t]j,n) and /^(x, < ; j, n), respectively, 
are defined in Eqs. (4.13) and (4.18). /^(x,f;j,n) and ft^(x,t ;ji,n) will be defined im- 
mediately. 

Both f m and / m ,jt are functions of tx m , m = 1,2,3. Let (/ m )" and (/ m ,J t)", respec- 
tively, denote the values of f m and f m ,k when u m , m = 1,2,3, respectively, assume the 
values of (u m )*, m = 1,2,3. Let 


(U; = £(/m,*)>(“**)" > 


m = 1 , 2 , 3 , 


fc=i 


and 


(/mt)? = £(/m, *)"(«*)", "» = W- 

fc=l 

Using an argument similar to that leading to Eq. (4.18), we assume that 

fc(x,t;j, n ) = (/»)? + (/.,);(* -ij) + (/-.)*(*-«■). "> = 1.2.3- 

As a result of Eq. (5.5), we also assume that 

9K(x,t;j,n) 


(5.9) 


(5.10) 


(5.11) 


hn(x,i;j,n)~ fe(x,t;j,n) - 


dx 


-, u^(x,t;j,n) , m = 1,2,3. (5.12) 


Also, we assume that, for any (x,t) £ SE (j,n), u m = ;j, n), f m = f^(x,t ;j,n), 

and / m = f^(x,t;j,n) satisfy Eq. (5.4), i.e., 


du^x,*; j,n) _5_ 
dt dx 


**, , • ^ dfm( x ,t ;j, n) 

f m ( x . * ; * n ) ^ 


= 0. 


(5.13) 


The above condition again leads to Eq. (4.21). Thus, the diffusion term in Eq. (5.4) has no 
impact on how u* m (x,t ; j, n) varies with time within SE (j, n). This same fact was observed 
in Sec. 2. The reason behind it and its significance were also discussed there. As a result 
of Eq.(4.21), and other definitions given earlier in this section, one can conclude that the 
only independent discrete variables needed to be solved in the current solver, as in the 
Euler solver described in Sec. 4, are also (u m )" and (u mi )”. 

A comparison between Eqs. (4.20) and (5.13) reveals that, for the current solver, 
Eqs. (4.22) and (4.23) should be replaced by 


&t>m(x,t;j,n) 

dt 




d/mOM;i,n) 

dx 


(5.14) 
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and 


(5.15) 


cfy m (r,t;j,n) 

dx 




respectively. Note that Eqs. (5.15) and (4.23) are identical. According to Eq. (5.11), the 
second term on the right side of Eq. (5.14) is simply the constant — (/mi)"- Thus, for the 
current solver, Eq. (4.24) should be replaced by 


tMM ; j,n) = (/„)?(* - n - (u m )?0r - xj) + (1/2 )(/„,)"(* - <") 2 
- (l/2)(u mt )"(x - X ,) 2 + (/„)?(* - xy)(( - f), 


(5.16) 


def 


where 

(/■)? = (/«)? - (/-.)"• (5-17) 

The only difference between Eqs. (4.24) and (5.16) is that ( / m )” in Eq. (4.24) is replaced 

* 

by (fm)j in Eq. (5.16). Obviously, Eq. (4.25) is still valid for the current solver. Because 

;j,n) is independent of (/ m )” and (/ m t)", Eq. (4.25) implies that the last two 
parameters are irrelevant in flux evaluation. Moreover, because the current solver will be 
constructed using only flux-balance conditions, these parameters are also irrelevant in the 
following construction. 

For all (j, n) £ ft, we assume that 


/ 

Jst 


h* m ■ ds = 0. 


’S{CE±(j,n)) 

With the aid of Eqs. (5. 16) and (4.25), Eq. (5.18) implies that, for all (j,n) £ ft, 

(«»)7 - (»»)£# ± " [(<*««)£$ + (»-.)?] 


(5.18) 


At 


AX 


[(/m)j ± ,/ 2 (/m)> ± \ [(^ m( )j±l/2 + (/"»<)>] — °- 


(5.19) 


Adding the two equations given in Eq. (5.19) results in 

1 r ' \"-l/2 , / 1/2 I / ' \ n — 1/2 


(llm)” — 2 [( u m)j_l/2 + ( U ”i)j + l/2 + (®m)j-l/2 (^m)j +1 / 2 ] 

where, for all (j, n) £ ft, 


/ ' \n def )n i ^ (t ^n , (^0 / r \n 

~ 4 (“*«*)> + AX Wm)j + (fmt)j, 


m = 1, 2, 3. 


(5.20) 


(5.21) 


Eqs. (5.20) and (5.21) are the current counterparts of Eqs. (4.28) and (4.29), respectively. 
By using Eq. (5.20), (u m )" can be solved explicitly in terms of discrete variables at the 
next lower time level. 
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By substraction of the two equations given in Eq. (5.19), and using Eq. (5.17), one 


has 


m = 1,2, 3, 


where, for all (j, n) G Q, and m = 1,2, 3, 


(M? = + § [MT+vl - 


1/2 


fi V * -1/2 


/" ')«- 1 / 2 


(5.22) 


(5.23) 


Note that (/ m )", m = 1,2,3, are functions of (u m )", m = 1,2,3, and the latter can be 
evaluated by using Eq. (5.21). Thus, ( 6 m )" , m = 1,2,3, can also be evaluated in terms of 
the variables at the (n — l/2)-th time level. 

Moreover, for all (j, n) G fl, let 


/ f+ \n 4«f 


4a* 

(ax ) 2 



m, k = 1, 2, 3, 


(5.24) 


and 


(««*)? = + (/;,*)? - E m >* = w - 


(5.25) 


i=i 


Then, with the aid of Eqs. (4.38)-(4.40) and (5.9), Eq. (5.22) can be reexpressed as 


3 


E (<*■»*)?«)? = < m ”. 

fc=i 


m = 1, 2, 3. 


(5.26) 


Because (/+ fc )” and (/**)", m, fc = 1,2, 3, are all functions of (u m )”, m = 1,2,3, so are 
(amt)? , rn,k — 1,2, 3. Thus, (a mJ t)” can also be evaluated in terms of the variables at the 
(n — l/2)-th time level. It follows that, for each (j,n) G fl, Eq. (5.26) represents a system 
of three linear equations for three unknowns (u^ x )”, m = 1,2,3. These unknowns (and 
thus (u mx )j, m = 1,2,3, through Eq. (4.39)) can be solved easily by a matrix inversion. 
Eqs. (5.20) and (5.26) form the current marching scheme. 

For all ( j , n) G fi, let (F + )” denote the matrix formed by (/^ fc )”, m, k = 1, 2, 3. Also, 
let 

(C±)” d = /±(F+)?, (5.27) 

and 

(At)” =' / - [(F + )"f ± (F + )", (5.28) 

where (F + )” is defined in Sec. 4. Note that (£>+)" is the matrix formed by (a m *)”, 

m,k — 1,2,3. Existence of its inverse [(Z) + )”] 1 will be assumed. As a resut, one can 
define 
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r\( \ def 1 

Q+U> n ) = 2 


[(£>+)?] -1 « 7 -)?(c + 


(.D-TjZl'il 


\ 


(5.29) 

-&-)%% \ 


and 


/ 


. def 1 

Q-0> n ) = g 




V [(B +);]' 1 (c + )j(C-)"-7 2 2 - [(X>+)J ] _1 (c + )"(B -)"-‘/ 2 2 } 


(5.30) 

Using Eqs. (5.27)-(5.30), and mathematical manipulations similar to those leading to 
Eq. (4.48), Eq. (5.20) and (5.26) can be cast into the matrix form 


q(j,n) = Q+ (j,n)q(j - 1/2, n - 1/2) + Q_ (j,n)q(j + 1/2, n - 1/2), 


(5.31) 


where q (j,n) is defined in Eq. (4.45). Note that q(j,n) is converted into q(j,n) (see 
Eq. (2.15)) if u” and (u+)y, are replaced by uj and (ax/4 )(u x )", respectively. Also 
Q+0» and Q_(j, n) axe converted into Q+ and Q-, respectively, if (i) (i r+ )" and 

(F+) 

^±1/2 all replaced by i /, and (ii) (F + )” and {F + ) n - ± l^ \ are till replaced by f. 
Thus, Eq. (5.31) is converted into Eq. (2.14) after the above substitutions. 


34 


l 



6. Stability Analysis 


The stability of the a-/z scheme will be studied using the von Neumann analysis. For 
all (j, n ) € ft, let 

q(j,n) = q k (n 1 0)e ij6 (i = v^, -tc < 0 < tt) (6.1) 

where g*(n, 0) is a 2 x 1 column matrix. Substituting Eq. (6.1) into Eq. (2.18), one obtains 

?(n + 1,0) = [Q(v, £, 0)] 2 9>, *) ( 6 . 2 ) 


where 


Q(v,Z,6) d = e~ i9/2 Q+ + e i9 > 2 Q-. 


(6.3) 

According to Eq. (6.2), the amplification matrix is the square of the matrix Q{v,(i,0). 
Substituting Eqs. (2.16) and (2.17) into Eq. (6.3), one has 


'cos( 0 / 2 ) — ii/sin( 0 / 2 ) 
QKf, 0 )=| i(l-u 2 )sm(0/2) 

i - v 2 + i 


-i'(l - V 2 - 0 sin(0/2) ^ 

1 ~ ^2 ~ ^ -[ cos (^/ 2 ) + ^sin( 0 / 2 )] j 


Let 


77 ( 1 /, £, 0 ) d = £cos( 0 / 2 ) -iv(l - i/ 2 )sin( 0 / 2 ). 
Then the eigenvalues of Q(u, £,6) are 


<T±(V,Z,0) = 


def £, 0 ) ± £, 0)P + (1 - t/2 ) 2 


i - v 2 + i 


Thus the amplification factors and G^_ ’ of the a-fi scheme are given by 

G ( ± = k±(^^^)] 2 - 


(6.4) 


(6.5) 


(6.6) 


(6.7) 


Note that 2 

G™ -> 1 and G ( l ] ( y~ 2 t| ) 33 6 0 (6 ‘ 8) 

if 1 - v 2 > 0. Because the amplification factor of a plane- wave solution to Eq. (2.1) 
approaches 1 as 0 -» 0 , G+* and G ( J } are referred to as the principal and the spurious 
amplification factors, respectively. Moreover, Eqs. (6.5)-(6.7) imply that 


G^ = 


|cos( 0 / 2 ) - xv sin( 0 / 2 ) ± yj[£ cos( 0 / 2 ) - it/sin( 0 / 2 )p 


-,2 


+ 1-^ 2 


i + e 


(6.9) 
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if 1 - i/ 2 ± 0, and f d = f/(l - i/ 2 ). Similarity between Eqs. (6.9) and (2.21) was noted in 
Sec. 2. 

In [1], the stability of the a-ji scheme is studied using a rigorous discrete Fourier 
analysis. The von Neumann stability analysis can be considered as a limiting case of the 
discreate Fourier analysis. By using Eqs. (4.33) and (4.34) in [1], one can infer that the 
a-\i scheme is stable if and only if, for all 8 with — 7 r < 8 < tt, 


max{|G+^|, |G^|) < 1 if Q(i/,£,0) is nondefective (6.10) 


and 

|C?^| < 1 if Q(v , £, 8) is defective. (6.11) 


Note that G+\ — G ^ if Q(v, £, 8) is defective [23, p.353]. Assuming £ > 0 and 1 — i/ 2 +£ ^ 0 
(the latter is a basic assumption of Eq. (2.14)), it is proved in [1] that the current scheme 
is stable if and only if v 2 < 1. 

Let (1 — i/ 2 ) 2 ^ £ 2 such that both Eqs. (2.14) and (2.24) axe valid. Combining 
Eqs. (6.5)-(6.7), one has 


G^G^ 



( 6 . 12 ) 


Because the amplification factors of the backward-marching scheme axe (G^) -1 and 

(G[? ) )~ 1 , stability of both Eqs. (2.14) and (2.24) requires that [G^ ) | = IGL^I = 1. Ac- 
cording to Eq. (6.12), the last condition cannot be met if fx > 0 and v 2 ^ 1. This result 
was used in a discussion given in Sec. 2. 

Next we study the stability of the a-e scheme. By substituting Eq. (6.1) into Eq. (3.7), 
one has 

9> + M) = [M(e, u, 0)] 2 g*(n,0) (6.13) 

where 

M(e, u, 8) d = e -,e/2 M+ + e' 0/2 M-. (6.14) 

According to Eq. (6.13), the amplification matrix of the a-e scheme is the square of the 
matrix M(e, v,8). Substituting Eqs. (3.8) and (3.9) into Eq. (6.14), one has 


M(e,M) = 


cos(0/2) — ius,\n(8/2 ) — i(l — i/ 2 )sin(0/2) 

i(l — e) sin(0/2) (2e — 1) cos(0/2) — ii/sin(0/2) 


(6.15) 


The eigenvalues A±(e, u, 8) of M(e, u, 8) were given in Eq. (3.13). The principal amplifica- 
tion factor G+* and the spurious amplification factor G^ ) of the a-e scheme were given in 
Eq. (3.12). Note that 


G+^ — ♦ 1 and G^ — > 2e — 1 as 8 —* 0 


(6.16) 
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if Eq.(3.14) is assumed. Moreover, from Eqs. (6.10) and (6.11), one infers that the a-e 
scheme is stable if and only if, for all 0 with — -k < 0 < ir, 

max{|G+ ) |, |G^|} < 1 if M(e, v,d) is nondefective (6-17) 

and 

|G+^| < 1 if M(e, u, 0) is defective. (6.18) 

Eq. (3.13) implies that 

|A + (e,^0)||A_(e,^0)| = |2e-l|. (6.19) 

By using Eqs. (3.12) and (6.17)-(6.19), one concludes that stability requires that |2e- 1| < 
1, i.e., 0 < e < 1. This is the first part of Eq. (3.14). Eq. (3.13) also implies that 

A±(e, v, 7 r) = —iv ± \/(l — e)(l — i/ 2 ). (6.20) 

Thus, 

max{|A + (e, y, 7r)|, |A_(e,v, 7r)|} > 1 if v 2 > 1 and e < 1. (6-21) 

The first part of Eq. (3.14) coupled with Eqs. (6.17), (6.18), and (6.21) implies that i/ 2 < 1 
is necessary for stability. Because the case t/ 2 = 1 is ruled out by the basic assumption 
1 - i/ 2 ^ 0 of Eq. (3.6), the second part of Eq. (3.14) is now proved. 

To prove Eqs. (3.17) and (3.18), note that Eq. (3.16) implies that 




X±{e,v, 0) = e (x' T x") 

(6.22) 

where 


x' d = (1 - i/ 2 )sin 2 (0/2) + 2(1 - e)cos 2 (0/2). 

(6.23) 

and 

x" 

d = 2cos(0/2)^/(l - e) [(1 - e)cos 2 (0/2) + (1 - i/ 2 )sin 2 (0/2)]. 

(6.24) 

With the aid of 

Eq. (3.14) and —ir < 0 < 7r, Eqs. (6.23) and (6.24) imply that 




x' = x" — 0 if e = 1 and 0 = 0, 

(6.25) 



, J = 0, if t = 1 and 0 = 0; 
X \ > 0, if e ^ 1 or 0 / 0, 

(6.26) 



X " > 2(1 - e)cos 2 (0/2) > 0, 

(6.27) 



x' ~ X" < (1 - ^ 2 )sin 2 (0/2), 

(6.28) 
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and 

(x' - x")(x' + x") = (xT - (x'T = (1 - i' 2 ) 2 sin 4 (0/2). (6.29) 

For the case e = 1 and 0 = 0, Eqs. (3.17) and (3.18) follow immediately from Eqs. (6.22) 
and (6.25). Thus, in the following proof of Eqs. (3.17) and (3.18), we assume that 

e / 1 or 0/0. (6.30) 

Combining Eqs. (6.26), (6.27), and (6.30), one concludes that 

x' + x" >0. (6.31) 

Eqs. (6.29) and (6.31) imply that 

x'-x">0. (6.32) 

Eq. (3.18) now follows from Eqs. (3.14), (6.22), (6.28), and (6.32). The validity of the 
first inequality sign in Eq. (3.17) follows from Eq. (3.18) and the fact that e(l — e) > 0 if 
0 < e < 1. The validity of the second inequality sign follows from the fact that 

X-(e, u, 0) - x+(e, v, 0) = 2ex" > 4e(l - e)cos 2 (0/2). (6.33) 

Eq. (6.33) is a simple result of Eqs. (6.22) and (6.27). To establish the validity of the last 
inequality sign in Eq. (3.17), note that 

X_(e, v, 0) = e(x' + x") = «[2x # - (x' - x")] < 2e*' 

= 2e [(1 - t/ 2 )sin 2 (0/2) + 2(1 - e)cos 2 (0/2)] , (6.34) 

< max{2e(l — v 2 ),4e(l — e)} < 4e 

where Eqs. (6.22), (6.32), (6.23), and (3.14) have been used. Moreover, because |G^| > 0, 
Eq. (3.16) implies that 

X-(e,^,0)<l. (6.35) 

The validity of the last inequality sign in Eq. (3.17) now follows from Eqs. (6.34) and 
(6.35). Q.E.D. 

Next we shall prove that Eq. (3.14) is also sufficient for stability. Note that, as a 
result of Eqs. (3.17) and (3.18), 0 < x±( e ? t/ )^)> ^d thus |G^| < 1, for all e, z/, and 0 
satisfying Eq. (3.14) and —it < 0 < 7r. As a result, Eq. (6.17) is always satisfied. To 
complete the proof, we need only to show that Eq.(6.18) is also satsfied. To proceed, note 
that = G^} if is defective. From Eqs. (3.12)-(3.14), one also concludes 

that e = 1 is necessary if = G_ . Moreover, Eq. (6.15) implies that M( 1, v, 0) is the 
identity matrix. Thus, one concludes that e = 1 and 0/0 are necessary if M(e, v, 0) is 
defective. Because (i) 

G±^ = [cos(0/2) — zVsin(0/2)] 2 if t = 1, (6.36) 
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and (ii) 

|[cos(0/2) — it/sin(0/2)] 2 | <1 if u 2 < 1 and 9 ^ 0, (6.37) 

one arrives at the conclusion that Eq. (6.18) is also satisfied. Q.E.D. 

Next we shall study the stability of the Euler solver defined by Eqs. (4.43) and (4.44). 
Because the von Neumann stability analysis, strictly speaking, is applicable only to a 
system of constant-coeeficient linear equations, we begin with a linearization of the above 
nonlinear equations. 

Recall that the independent variables for Eqs. (4.43) and (4.44) are (u m )" and (u ^ nx )^ , 
with m = 1,2,3 and (j,n) 6 ft. The other variables in these equations are functions of 
them. Let (u m )" and (u+ r )" also denote a solution to these equations. Let (u m )” + ^(«m)" 
and (u+ r )y 4- £(«m X )" denote another solution with 6(u m )” and <5( u mr)/ being small 
perturbations to (u m )" and (u+ r )", respectively. One may consider the second solution 
(hereafter to be referred to as the perturbed solution) to be the result of the first solution 
(hereafter to be referred to as the background solution) being perturbed initially by round- 
off errors at some time level. The purpose of the stability study is to determine whether 
the induced perturbation will amplify or die off as it propagates down the subsequent time 
levels. 

In the current linearization, we assume that 


|(«m*)>| < |(“ m )"|, for m = 1,2,3, and (j,n) £ ft. (6.38) 

According to Eqs. (4.13) and (4.39), the above assumption is equivalent to 

K,(xj + ax/ 4, t n ;j,n) - u^(xj,t n < \u* m (x j ,t n ;j,n)|, (6.39) 

i.e., the change in the value of u* m in SE(j, n) over a spatial distance of ax/4 is negligible 
compared with the value of at (xj,< n ). Similarly, we assume that 


, \n — 1 /2 / \n — 1/2 I „ 1/ \n-l/2 

( u m)j + 1 / 2 ( u m)j_ i /2 | ^ |( Um )j±l/2 

With the aid of Eq. (4.33), Eq. (6.40) implies that 


(«-)£$- (a-)? « !(*«)"! 


(6.40) 


(6.41) 


Because the magnitude of the round-off error of a small quantity is not necessarily 
smaller than that of a large quantity, in contrast to Eq. (6.38), we assume that 

K( u mr)j | ~ |<5( u m)"| , for m = 1,2,3, and (j,n) e ft. (6.42) 

Here the symbol implies that the quantities on both sides have the same order of 
magnitude. Also because round-off errors could vary erratically from one mesh point to 
another, in contrast to Eq. (6.40), we assume that 


^m) 


n — 1/2 
J + l/2 


6(u m ) 


n-1/2 
3 - 1/2 



n-1/2 

i±l/2 



(6.43) 
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Let 


denote the perturbed- solution counterpart of (/+ * )>±l/2( Ufc )>±l/2' Then 

t 6 - 44 ) 


The proof follows directly from the fact that , m, A: = 1, 2, 3, are homogeneous functions 
of degree 0 [20, pll] in the variables u m , m = 1,2, 3, and thus 


Y' 9 2 fm 

" du k dui Ul 


= 0, 


m,k = 1,2,3. 


(6.45). 


Q.E.D. With the aid of Eqs. (6.38) and (6.40)-(6.44), linearization of Eqs.(4.43) and (4.44) 
can now proceed by using the fact that both background and perturbed solutions satisfy 
these two equations. 

Recall that, as defined in Eqs. (4.7) and (4.8), m, k = 1, 2, 3, are functions of u m , 
m = 1,2,3. In Sec. 4, we also define (/ m> *)" to be the value of / m , * when u m , m = 1,2,3, 
respectively, assume the values of (u m )", m = 1,2,3. Moreover, (u m )", which is different 
from (u m )", is also defined in Eq. (4.33). To simplify the linearized versions of Eqs. (4.49) 

and (4.44), we define (/ m ,Jt)> to be the value of / m> * when u m , m — 1,2,3, respectively, 
assume the values of (u m )", m = 1,2,3. Furthermore, let 


( i+ \n de* 


m, k = 1,2,3. 


(6.46) 


Then the linearized versions of Eqs. (4.49) and (4.44) can be expressed as 


^ El «-* + (/£*>”! 


k~i 


+ 




i=l 




+ - (/+*)?] sk)”;, 1 # 


- E(/:.,)"(/A)" 


/=1 


r/ 4- 

^( U *r)j+l/2 


}■ 


(6.47) 


m = 1, 2, 3, 
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and 


«(•£,)? = \ e { [(*-»)? - <«*] [<(<*»);:,’ 'il - «(»*);+$] 

Z *=1 l 

+ [2(e mt ); - d mi + (/+,*)?] S(ui,)'r}% 

+ [2(e„oy - 6 m „ - (/+,,)?] «(<&)£$ }. 


m = 1,2, 3, 


(6.48) 


respectively. Note that, in arriving at the final forms given above, we have replaced both 
(/£*)“+ $ and with (/+*)?. According to Eq. (6.41), these substitutions 

introduce only errors which are higher order than every terms presented in Eqs. (6.47) and 
(6.48). 

A comparison among Eqs. (4.43), (4.44), (6.47), and (6.48) reveals that the first two 
equations will turn into the last two equations if, for all m, k = 1,2,3, and O', n) £ ft, 
(«•»)?, («+,)? and (/+*)£$ are replaced by «5(u m )", «5(u+ r )" and (/+*);, respectively. 
Roughly speaking, one can say that Eqs. (6.47) and (6.48) can be obtained from Eqs. (4.43) 
and (4.44) by “freezing” the coefficients and replacing the “background” variables (u m )” 
and (u* r )“ with the “perturbation” variables S(u m )” and S(u ^ lx )^ , respectively, for all 
m = 1,2,3, and O', n) £ ft. 

The von Neumann stability analysis for Eqs. (6.47) and (6.48) can be performed easily 
after each of them is decoupled. To proceed, first we convert them into matrix forms. For 
all 0, n) £ ft, let 6u” and <5 (tT+)", respectively, denote the column matrices formed by 

S(u m )] and S(u+ X )], m = 1,2,3. Let (F+)?, (G)", and (G _1 )"> respectively, denote the 
3x3 matrices formed by (/**)", (dmk)", arid ( 9m \ )?’ rn,k = 1,2,3. Then Eq. (4.34) is 
equivalent to 

(£)? = (G)7 Hag ((«,)?,(« 2 )7, («.);) (<?-■)”. (6.49) 


where (£)" was defined in Sec. 4. Furthermore, Eqs. (6.47) and (6.48) can be reexpressed 
as 



I + (F + )J 


,--. 71 - 1/2 . 


+ 


7-(F + )" 


c -*n — 1/2 

^i+l/2 


[/-((F+)”) 2 ] «(«+);-■£ 

[/-((f + ) j “) 2 ]«(sJ)"; 1 1 / /2 |, 


(6.50) 
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and 


*( 


r-*n — 1/2 — 1/2 

S-l/2 “ S + l/2 


+ [2(B)? - I + (F + )”] «(«£)£$ + [2(E)] - I - (F*)j] 6(Si )'-% } 

respectively. 

A result of Eqs. (4.12) and (4.40) is 

(G-')](F + )] (<?)“ = dia g ((P t )],(v 2 )], (»„)]), 

where 

(e , 9 if!2z3ki, ^ 2£. 

AX AX AX 

By using Eqs. (6.49) and (6.52) and matrix manupulations such as 

= (g-')](f+)](6)](g-')]6u]-;i 2 2 , 

and 

2 1 /<“ll / A . A. A \ 2 A 


(6.51) 


(6.52) 


(6.53) 


(6.54) 


(G” 1 )? ((*+)*) *(«+)£# = (((T^F+^G)?) (G-'W&VZlil (6-55) 


it is easy to see that both Eqs. (6.50) and (6.51) can be decoupled if the expressions on 
both sides of them axe multiplied from left by the matrix (G -1 )". For all (j, n) £ ft, let 
6'(u m )j and S'(u^ x )J, m = 1,2,3, respectively, denote the components of the column 
matrices (G -1 )” <5u” and (G -1 )" 8{u % )”. Then the decoupled equations can be written as 


«'(“•»)" = j{ [i + (*»)?] 8 + [i - ((->™)?) 2 ] tt’te,)':}* 

+ [1 - (*„)"] S'(u m )l;$ - [l - «i> m )") 2 ] «■(«+, )£$}, m = 1,2,3, 


(6.56) 




+ (2(e„); - 1 + (£>„);] 6'(uL)]:% 

+ [2(l m )J - 1 - (i>„ )?] «>+,)£$}. m = 1,2,3. 


(6.57) 
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With the aid of Eqs. (2.15), (3.8) and (3.9), one can see that, for each m, Eqs. (6.56) 
and (6.57), respectively, will be converted into the two component equations contained 
in Eqs. (3.6) if, for all (j, n) G 12, 6'(v m )*, 6'(u+ x )?, (t> m )>, and («m)", respectively, are 
replaced by u”, (ax/4 )(«*)", v, and e. Moreover, if we consider only the Fourier compo- 
nents of 6'(u m )j and <5'(u+ r )” with sufficiently short wavelengths, i.e., and (e m )" do 

not vary substantially over these wavelengths, the last two coefficients can be considered 
as constants in the von Neumann stability analysis. Note that round-off errors generally 
sure dominated by short-wavelength Fourier components. As a result of the above consid- 
erations, at least approximately, one can obtain the stability conditions of the marching 
scheme defined by Eqs. (4.43) and (4.44) by a straightforward generalization of Eq. (3.14), 
i.e., the marching scheme is stable if, for all (j, n) G 12, 

0<(e„)”<l, and [(•>„)?] 2 < 1, m = 1,2,3. (6.58) 

By using Eq. (6.53), it is easy to see that Eq. (4.49) is equivalent to Eq. (6.58). 
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7. Consistency and the Truncation Error 


Consistency and the truncation error of the a-jx scheme were studied and given in 
Sec. 6 of [1]. In this section, a similar study for the a-e scheme is presented. 

As a preliminary, note that Eq. (3.7) can be expressed explicitly as 

u i +1 = + V ^ e + ( 2 “ + (V2)(l - v 2 )(e + v )± *(«,)?_! 

+ 2(2 - e)(l - u 2 )u] - i/(l - u 2 )ax(u x )] (7.1) 

+ (1 - u)[e - (2 - e ) v ] u ] +1 - ( 1 / 2)(1 - i/ 2 )(e - u)ax(u x )] +1 }, 

and 

A X(U X )] +1 = 

2(e - l)(e + i/)u"_ x + (1/4) [e(4e - 3) + 2(2e - l)u + (2 - e)i/ 2 ] a x(u x ) 1 -_ 1 

(7.2) 

+ 4(1 - e)t/it" + (1/2) [4e 2 - 5e + 2 + (e - 2)u 2 ] ax(u x )" 

+ 2(1 - e)(e - u)u] +1 + (1/4) [e(4e - 3) - 2(2c - l)v + (2 - e)u 2 ] A x(u x )] +1 . 


Eqs. (7.1) and (7.2) represent a system of two discrete equations for each (j, n) 6 ft. It will 
be shown in this section that a solution to a pair of particular partial differential equations 
(PDE’s) will satisfy the above discrete equations under certain limiting conditions. One 
of these PDE’s is Eq. (2.22). 

To proceed, let u(x, t ) and u(:r, t ) be two smooth functions, let 


w(x,t) = V (x,t) — 


du(x,t) 

dx 


(7.3) 


For all ( j,n ) € ft, let u" d = u(xj,f"), v " d = v(x j,t n ), and w” = f w(xj,t n ). Let 


it- 1 

[DEI]" d = -±j- - — 1(1 + u)[e + (2 - e)u]u ]_ x + (1/2)(1 - u 2 )(e + 

+ 2(2 - c)(l - v 2 )^ - i/(l - v 2 )axv ? 

+ (1 - v){e - (2 - e)i/]u? +1 - (1/2)(1 - v 2 )(e - u)axvJ +1 }, 


(7.4) 



and 


[DE2]" =' «; +1 

- (2/ax)(c - l)(e + - (1/4) [e(4€ - 3) + 2(2e - l)x + (2 - e)^ 2 ] 

(7.5) 

- (4/ax)(1 - e)uu] - (1/2) [4e 2 - 5e + 2 + (e - 2)v 2 ] v] 

- (2/ax)(1 - e)(e - i/)fi; +1 - (1/4) [e(4e - 3) - 2(2e - l)i/ + (2 - e)v 2 ] v? +l . 

Note that [DEI]™ and [DE2]] axe defined such that Eqs. (7.1) and (7.2) axe equivalent to 

[DEI]] = 0, (7.6) 

and 

[DE2]] = 0, (7.7) 

respectively, if, for all (j, n) € fl, u” and v] in Eqs. (7.6) and (7.7) are replaced by u” and 
(u x )", respectively. 

Substituting Taylor series expansions of u" +1 , Uj±i, u” +1 and u" ±1 about (xj,t n ) into 
Eqs. (7.4) and (7.5), one concludes that, for any smooth functions u(x,t ) and i/(x,t), 

[DEI]] - [J PDE ]] = [ER1]], (7.8) 


[DE2]] - 4e(l - e)w] = [EB2]]. 

Here, assuming all derivatives are evaluated with x = xj and t = t n , 

, , <9£t 5u 

+. 


ndet Atfd 2 U 2 d 2 u\ f f(Ai ) 2 


^ J — 2 \ dt 2 Q dx 2 / 4 [ a< J dx 

(At ) 2 /d 3 u 3 d 3 u\ a [(ax ) 2 - a 2 (At) 2 ] /^u d 2 u> 

6 V Ht 3 + Q dx 3 / 24 [dx 3 dx 2 


(At ) 2 

/ d 3 u 

6 

V dt 3 

(At ) 3 

(d*u 

24 

V dt 4 

1 

(AX ) 2 

48 [ 

At 


i\ (ax ) 3 /ax \ d 3 w 

0 + ^r € \~zi~ a )d^ 


is ■ 2a2 H [(Ax)2 _ q2(a<)2 } k + °[ (ax) 1 + o[At(Ax) 3 ] 

0[(a*) 2 (ax) 2 ] +0[(At) 4 ] +e|o[At(Ax) 3 ] + ^0 [(ax) 4 ]}, 


(7.10) 
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and 


- i [e(4e - 3 )(ax) 2 + (2 - «)a 2 (Af) 2 ] 0 

d / d 2 u 

di V 


+ 


( A *) 2 

2 

(A if 


d 2 w 

~W 

d^w 

~W 


+ 


dt 2 ° dx 2 


d_ (cPu 3 guV 
dx \ dt 3 ' dx 3 ) 


dw 

dx 


e [(4c - 1)(ax) 2 - 3g 2 (a<) 2 ] d z u 


12 


dx 3 


+ 


(2e — l)aAf(Ax) 2 cPw 
6 dx 3 


+ 

+ c 


aAt [c(ax) 2 — a 2 (A<) 2 ] d 4 u 

6 


dx 4 +0[a<(ax) 3 ] +0[(a*) 2 (ax) 2 ] + 0[(a<) 4 ] 
|o[(ax) 4 ] +0[a^(ax) 3 ] + 0[(a<) 2 (ax) 2 ] | + c 2 0[(ax) 4 ] . 


(7.12) 


Because we assume that the parameter e can vary with At and ax, it is not treated as a 
constant in the derivation of Eqs. (7.8)-(7.12). In other words, a 11 the order-of-magnitude 
quantities given in Eqs. (7.11 ) and (7.12) axe independent of e. 

The significance of Eqs. (7.8) and (7.9) will be discussed under different assumptions 
about e. Assuming that c(l — e) ^ 0, Eq. (7.9) can be rewritten as 



[DE2']" - w’j = [ER2']j, 

(7.13) 

where 



[DE2 1 ]” d = [DE2] 

7/|4<l - e)], and |£K2’|" d = [M2]"/[4 e (l - <)]. 

(7.14) 

The following comments 

are made for Eqs. (7.8) and (7.13): 


(a) For all ( j, n) € ft, 

[PDE]j = 0, and w n } = 0, 

(7.15) 

if u = u(x , t ) and v 

= v(x,t) satisfy both Eq. (2.22) and 



du 

v- — = 0. 

ox 

(7.16) 


(b) Note that Uj = and (u x )" = v" satisfy Eqs. (7.1) and (7.2) if u r - and vj satisfy 
Eqs. (7.6) and (7.7). As a result, [DE 1]" and [DE2'}^ can be considered as the a-e 
scheme’s approximations for [PDE}J and u>", respectively. Eqs. (7.8) and (7.13) then 
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state that [EDI]” and [ED2']" are the errors of these approximations, respectively. Let 
[EDI]? -+ 0 and [ED2']? -» 0 as At, ax -+ 0. Then [DEI]] and [DE2']?, respectively, 
approach [ PDE ]] and w] as At, ax — * 0. Note that the limits [P DE]] and w] are 
independent of At and ax. 

(c) With the aid of the observations made in (a) and (b), one concludes that Eqs. (7.1) 
and (7.2) may be considered as the discrete approximations of Eqs. (2.22) and (7.16), 
respectively, with the understanding that u ] and (u z )” are the discrete counterparts 
of u and v, respectively. Note that Eq. (7.16), i.e., v = du/dx, is consistent with the 
fact that (u x )j is the numerical analogue of the value of du/dx at the mesh point 

(Xj,t n ). 

(d) Let u = u(x,t) and v — v{x,t ) be a solution taEqs. (2.22) and (7.16). Then Eq. (7.3) 
implies that w(x,t) = 0. Moreover, it can be shown that 


d t u(x,t) ( ^ d e u(x, t) _ n 

dt t { ’ dx l 

As a result, Eqs. (7.11), (7.12) and (7.14) imply that 


£ = 1,2,3, 


(7.17) 


a [(ax) 2 - a 2 (At) 2 l cPu 

[£*i] j.JL-L—L-LL 


+ 
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( ai)2 _2a ! A t 


At 


[(ax) 2 - a 2 (a() 2 ] ^l+O^Axfl + OlAttAx) 8 ] (7.18) 


ax 4 


+ 0[(a<) 2 (ax) 2 ] +0[(a*) 4 ] + e{o[A<(Ax) 3 ] + ^-0 [(ax) 4 ]}, 


and 


rriiwlin [(4e — 1 )(ax) 2 — 3a 2 (At) 2 ] cPu t aAt [{(ax) 2 - a 2 (At) 2 ] d*u 
[ER2 ]] = 48(1 - e) + 24e(l - e) dx* 

+ — i-^{o[A((Ax) 3 ] +0[(a() 2 (ax) 2 ] +0[(a<) 4 ]} (7.19) 

+ J^ T ^y{0[(AX ) *] +0[a((ax) 3 ] +0[(A«) 2( AX) 2 ] +.0[ ( AX) 4 ]}. 

As in Eqs. (7.10)-(7.12), the derivatives in Eqs. (7.18) and (7.19) are evaluated with 
x = Xj and t = t n . 

(e) Again, let u = u(x,t) and v = v(x,t) be a solution to Eqs. (2.22) and (7.16). Then, 
by definition, [DEI]] and [DE2']] are the truncation errors of Eqs. (7.1) and (7.2) 
with respect to the above PDE’s [24, p.20]. Futhermore, we have [PDE]] = 0 and 
w] — 0. Thus 

[DEI]” = [EDI]", and [DE2']” = [ED2']" , (7.20) 
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i.e., [Ei?l]” said [ER2']j given in Eqs.(7.18) and (7.19) axe the truncation errors. 

(f) Consider the special case in which t does not vary with At and ax. According to 
Eqs. (7.18), and (7.19), we have 


[ER1]] - 


c(ax) 4 

48a< 


d*u 

dx 4 


+ 0(ax) 


0 and [ER2'}” -► 0, 


(7.21) 


as At, ax — > 0, regardless how the mesh is refined. Thus, [Ei21]” — ► 0 as At, ax — ► 0 
only if the mesh refinement is subjected to the condition: 

(ax) 4 /a< — ♦ 0 as At, ax — * 0. (7.22) 

Thus, the a-t scheme is consistent with Eqs. (2.22) and (7.16) if the mesh refinement 
is subjected to Eq. (7.22). 

Note that, by using Eqs. (7.8), (7.13) and (7.21), it can be shown that the a-t 
scheme is consistent with Eq. (7.16) and the modified equation [25-27]: 


du du 
dt + a dx ~*~ 


e(Ax) 4 d 4 U € r 51 


(7.23) 


if Eq. (7.22) is not satisfied. Eq. (7.23) differs from Eq. (2.22) in the presence of a 
leading diffusion term and other higher-order terms. As a result, for a constant t, the 
a-t scheme becomes more diffusive as the ratio (ax) 4 /a t increases. 

(g) In the general case in which t = e(At, ax) , At > 0 and ax > 0, consistency of the a-t 
scheme is dependent on the behavior of e(At,Ax), as At, ax — * 0. As an example, let 
«2 and ei be two constants such that 62 > t\ > 0. Let 


e 2 > |e/ At] > ti as At, ax — »• 0. (7-24) 

Then Eqs. (7.18) and (7.19) imply that the a-t scheme is consistent with Eqs. (2.22) 
and (7.16). As another example, let e<j > 0 be a constant, and 

f AT 

> t 0 as At / ax — ► 0. (7.25) 

At 

Then the a-t scheme is also consistent with Eqs. (2.22) and (7.16) if the mesh refine- 
ment is subjected to the condition 

At/ ax — ► 0 as At, ax — * 0. (7.26) 


This completes the discussion on Eqs. (7. 8) and (7.13). For either t = 0 or t = 1, the 
second term on the left side of Eq. (7.9) vanishes. For t = 0, which is a special case of the 
a-fi scheme with n — 0, it is shown in [1] how consistency of the a-t scheme can be studied 
by recasting Eq. (7.9) into another form. By applying the same technique, consistency of 
the scheme with t = 1 can also be studied. 
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8. Numerical Results 


In [1], numerical solutions of Eq. (2.1) generated by the MacCormack [4, p.102], the 
Leapfrog/DuFort-Frankel [4, p.161], and the a-// schemes are compared with the corre- 
sponding analytical solutions for different values of physical coefficients, mesh parameters 
and total running times. These comparisons show that the a-/i scheme is fax superior to 
the Leapfrog/DuFort-Frankel scheme in accuracy, and has a substantial advantage over 
the MacCormack scheme in both accuracy and stability. 

In this section, accuracy of both the Euler and the Navier-Stokes solvers will be 
evaluated numerically using a shock tube problem suggested by Sod [28]. Because the 
a-e scheme may be considered as a special case of the Euler solver, no separate numerical 
evaluation for the a-e scheme will be given. 

Let the specific heat ratio 7 = 1.4. At t = 0, let (i) (p,v,p) = (1, 0, 1), i.e., 
(itj , U2, U3) = (1, 0, 2.5), if x < 0, and (ii) (p,v,p) = (0.125, 0, 0.1), i.e., (ui,u 2 ,u 3 ) = 
(0.125, 0, 0.25), if x > 0. For all (j,n) 6 D, let xj = j ax, and t n = nAt. Then (i) 


(1, 0, 2.5), 


if j = -1/2, -3/2,...; 


;,(u 2 );,(u 3 ©= (8.1) 

[ (0.125, 0, 0.25), if j = 1/2, 3/2, . . ., 

and (ii) («„,*)$ = 0, j = ±1/2, ±3/2, . . ., for m = 1, 2, 3. Hereafter, we assume n > 0. 

The above initial conditions coupled with several equations given in Secs. 4 and 5, 
imply that, for both the Euler and the Navier-Stokes solvers, (u m )" is a constant and 
(«mr)y = 0 in two separate regions that are defined by j < — (n + 1/2) and j > (n + 1/2), 
respectively. Thus, one needs to evaluate the above variables only if (n -f 1/2) > |;|. 

Without exception, ax = 0.01 is assumed in this section. Also, all numerical results 



will be compared with the exact weak solution at t = 0.2. Because, at t — 0.2, the effect of 


the initial discontinuity at t = 0 is far from reaching the spatial regions defined by x > 0.5 
and x < -0.5, respectively, numerical computations will be simplified by assuming that, 
for all n with t n < 0.2, (i) 


f (1, 0,2.5), if xj < —0.5; 

((« 1 ©(« 2 ©(u 3 )7)= { (8.2) 

[(0.125,0,0.25), if Xj > 0.5, 

and (ii) (u m r)" = 0 if |x_,| > 0.5. Because ax = 0.01, the above assumptions imply that 
the computation domain can be limited to |j| < 50. 

In the initial evaluation, we consider the Euler marching scheme defined by Eqs. (4.28) 
and (4.32). numerical results (triangles) obtained assuming A t = 0.004 and e = 1/2 are 
compared with the exact solutions (solid lines) in Fig. 6. Because each marching step 
advances the solution from t to t + At/2, these results at t — 0.2 are obtained after 100 
steps. Also it can be estimated that CFL = 0.88 where CFL is defined to be the maximum 
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value of (|u| + |c|)a</ax. Thus the numerical calculation is carried out within the stability 
limits given by Eq. (4.49). Note that the agreements between the numerical results and the 
exact solutions axe excellent. Particularly, shock discontinuity is resolved almost within 
one mesh interval, and contact discontinuity is resolved in four mesh intervals. Also, there 
are only slight numerical overshoots and/or oscillations near these discontinuties. 

According to the discussions given in Secs. 3, 4, and 6, the Euler solver behaves like 
the Leapfrog scheme if c = 0, and like the Lax scheme if e = 1. The former is free from 
numerical diffusion while the latter is highly diffusive. The current scheme with e = 1/2 
can be considered as a scheme midway between the above two celebrated schemes. 

Moreover, the last term on the right side of Eq. (4.32) vanishes if e = 1/2. The 
remaining term is simply a central-difference approximation for (ti mi )”. 

Let Eq. (4.32) be modified with (u£, r )” being replaced by (uj£x)? (see eqs. (4.51) and 
4.54)). Again assuming that At = 0.004 and c = 1/2, the numerical results obtained with 
a = 1, or = 2, and a = 3, respectively, are given in Figs. 7-9. The effectiveness of the 
above modification as a tool to surpress numerical wiggles near discontinuities is apparent. 
It was explained in Sec. 4 why this modification does not cause the smearing of shock 
discontinuities. Furthermore, the modification has no discemable effect on the smooth 
part of the solution. Because (uj£° x )" = (u c mx )j if a = 0, in the following discussion, it 
should be understood that the above modification is turned off if a = 0. 

Note that the results shown in Figs. 6-9 can be generated using the sample program 
listed at the end of the present paper. It is coded assuming e = 0.5. The value of the input 
parameter ic is equal to that of a. 

Let a = 0 and At = 0.004. The numerical results obtained with e = 0.1, e — 0.3, 
e = 0.7, and e = 0.9, respectively, are given in Figs. 10-13. Note that the case with e = 0.5 
was presented in Fig. 6. For e = 0.1, because the scheme has very small numerical diffusion, 
pronounced wiggles appear in large regions near discontinuities. However, because of the 
same reason, the smooth part of the solution is highly accurate. For e = 0.3, the wiggles 
cure less pronounced and appear in more limited regions. Also the smooth part of the 
solution becomes less accurate. As the value of e increases, the wiggles disappear and the 
solution becomes more diffusive. The solution obtained with e = 0.7 is excellent except 
that, compared with the case with e = 0.5, it requires one more mesh interval to resolve 
the contact discontinuity. The results shown in Figs. 6 and 10-13 are consistent with the 
theoretical prediction that the Euler solver becomes progressively diffusive as the value of 
e increases from 0 to 1. 

Figs. 14 and 11 are both generated using the same conditions except that a = 2 in 
the former, while o = 0 in the latter. Note that the wiggles almost completely disappear 
in Fig. 14. 

The above numerical results Eire all generated assuming At = 0.004. The numerical 
results shown in Fig. 15 are generated with At = 0.002, e = 1/2, and a = 0. It is 
obtained after 200 steps and CFL = 0.44. A comparison between Figs. 6 and 15 reveals 
that the current solver is more diffusive at a smaller CFL. Note that, by considering the 
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truncation error, it was shown in Sec. 7 that, for constant e and ax, the a-e scheme 
becomes more diffusive as At decreases. A similar conclusion can also be reached by- 
studying the amplification factors given in Eqs. (3.12) and (3.13). Because the Euler 
solver is a straightforward extension of the a-e scheme, one would expect that the former 
also behaves similarly. Fig. 16 shares with Fig. 15 the same defining conditions except that 
a = 1 for the former. 

The numerical results shown in Figs. 17 and 18 are generated assuming A f = 0.0004 
(i.e., CFL = 0.088). Note that it takes 1000 marching steps to advance the solution to 
t = 0.2. Other defining conditions for these figures are identical to those for Figs. 6 and 
7, respectively. As expected, the results obtained with low CFL are much more diffusive 
than those obtained with CFL closer to 1 (see Figs. 6 and 7). Also, as the value of CFL 
decreases, the diffusive effect of replacing a = 0 with a = 1 becomes more discerable even 
in the smooth part of the solution. In other words, numerical diffusion introduced by 
replacing a = 0 with a > 0, is greater when CFL is small. 

To modify the above Euler solver such that it can compensate for the observed effect of 
increasing numerical diffusion as At decreases, in the following discussions, we shall consider 
the more general marching scheme defined by Eqs. (4.28) and (4.36). The parameter (e)” 
in Eq. (4.36) will be dependent on the mesh position (/, n ) and the ratio At/ ax. Moreover, 
the term (u^)" in Eq. (4.36) will be replaced by (u^)”, which is defined in Eq. (4.59). 
The weight factor j3 will also be dependent on (j, n) and At/ ax. 

To proceed, let 

£(x) *= f xexp(l — x), 0 < x < 1. (8-3) 

Because £ is an increasing function within its domain, we have 

C(x)<C(l) = l, 0 < x < 1. (8.4) 

For all (j, n) € 0, let 

(e)? = K(0w)”), (8-5) 

and 

«,)J = W (u»,+ )7;a, ) , (S.6) 

where (i) mai )" is defined in Eq. (4.50), and b and a are constants that do not vary from 
one mesh point to another. Because (e m )? = (e)”, m = 1,2,3, is assumed in Eq. (4.36), 
Eqs. (4.49), (8.4) and (8.5) imply that (i) (i> m ar)” is in the domain of £(x), and (ii) 
0 < b < 1. 

Note that (Umax)] is proportional to At/ ax. Thus, Eqs. (8.3) and (8.5) imply that (e)” 
is an increasing function of At/Ax, i.e., it decreases as At decreases if other parameters 
are held constant. Because numerical diffusion decreases as (e)” decreases, with other 
factors being equal, the replacement of a constant e with (e)” has an effect in reducing 
numerical diffusion a s At decreases. This effect will compensate for the observed opposite 
effect on numerical diffusion as a f decreases with e, ax, and the toted running time being 
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held constant. Furthermore, because £(x)/x — * 0 as x —* 0, Eqs. (4.50) and (8.5) imply 
that 

(l)?^^6(|C"l + |cJ|) as Ai/AX-^0. (8.7) 

Eq. (8.7) is similar to the consistent condition given in Eq. (7.25). 

Moreover, for a fixed a, W(x_, x+; a, 0) -+ (x_ + x + )/2 as 0 — ♦ 0. This fact coupled 
with Eq. (4.53) implies that the numerical diffusion introduced as a result of replacing 
( u mr)? with (u£ r )" will decrease as 0 decreases. Because (i> ma i)” is proportional to 
At /ax, with other factors being equal, the replacement of («„*)" by («£*)? defined in 
Eq. (8.6), has an effect in reducing numerical diffusion as At decreases. This effect will 
compensate for the observed opposite effect on numerical diffusion as At decreases with 
a, 0, ax , and the toted running time being held constant. Note that Wo(x_,x+;a) is a 
special case of W(x_, x + ; a, 0) with 0 = 1. 

Assuming a = 1 and b = 0.5, the numerical results shown in Figs. 19, 20 and 21 
are generated with At = 0.004 (CFL = 0.88), At = 0.0004 (CFL = 0.088), and At = 
0.0001 ( CFL = 0.022), respectively. Note that the results shown in Fig. 19 axe almost 
identical to those shown in Fig. 7 which were generated assuming the same conditions but 
using a simpler marching scheme. However, the results shown in Fig. 20 axe fax less diffusive 
than their counterparts shown in Fig. 18. One can conclude from this comparison and the 
results shown in Fig. 21 that the current modified Euler solver is capable of generating 
accurate numerical solutions even for the case with a very small CFL. 

In the above modified Euler scheme, (e)” and 0 are expressed as two special func- 
tions of (i>mar )” 5 respectively. They are only two among many possible choices. The 
investigation of other choices is a subject to be studied in the future. 

The most general marching scheme presented in Sec. 4 is that defined by Eqs. (4.28) 
and (4.35). It requires several matrix multiplications at each mesh points and, therefore, is 
much more costly. Thus, its use is difficult to justify unless a substantial gain in accuracy 
can be made. How this most general marching scheme can be applied wisely is left for a 
future study. 

This completes the numerical study of the Euler solver. We conclude this section with 
a numerical evaluation of the Navier-Stokes marching scheme defined by Eqs. (5.20) and 
(5.26). Again the initial conditions defined in Eq. (8.1) are assumed, and the numerical 
solutions are compared with the exact weak solution of the Euler equations at t = 0.2. 
The numerical results shown in Figs. 22-28 are generated assuming At = 0.004, ax = 0.01, 
7 = 1.4, and Pr = 0.72. The value of the Prandtl number used here is that for air at 
standard conditions. The values of the Rcl for these figures are 2,000, 4,000, 6,000, 8,000, 
10,000, 12,000, and 20,000, respectively. 

From the results shown in these figures, one concludes that, for a high-Reynolds- 
number flow, the shock can be resolved within one mesh interval by the current Navier- 
Stokes solver. Also the contact discontinuity can be resolved within a few mesh intervals. 
Note that these results are obtained without using any ad-hoc parameters or techniques. 
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Because the Reynold number is inversely proportional to the physical viscosity, as expected, 
numerical overshoots and oscillations shown in these figures increase slightly as the values 
of the Reynolds number increase. 

Furthermore, through repeated numerical experiments using different physical and 
mesh parameters, it is established that the current Euler solver is stable if, for all (j, n) € ft, 

0 < ReL, 0 < Pr, and (£ m0 r)? < 1 (8.8) 

However, because a Navier-Stokes problem is fundamentally an initial- value/boundary - 
value problem, the current explicit marching scheme obviously cannot model such a prob- 
lem unless the boundary effect is small, i.e., when the contribution of the viscous terms 
to Eqs. (5.20) and (5.26) is small compared to that of the convection terms. In general, 
this implies that the current scheme is applicable only to high- Reynolds-number flows. 
Note that the Leapfrog/Dufort-Frankel and the a-fj, schemes [1] also encounter a similar 
limitation in modelling Eq. (2.1). 

Finally, note that the current Navier-Stokes solver with ReL — oo (i.e., the physical 
viscosity vanishes) and Pr = 0 can be considered as a nonlinear extension of the inviscid 
a-fi scheme. Because the latter scheme is neutrally stable, generally one would expect that 
a nonlinear extension of such a scheme is unstable. However, it has been shown numerically 
that the current Navier-Stokes solver is stable even for the above limiting case as long as 
(j>mai)” < 1 for all (;', n) 6 ft. 
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9. Conclusions and Discussions 


Several key limitations of the finite difference, finite volume, finite element, and spec- 
tral methods were discussed in Sec 1. The method of space- time conservation element and 
solution element was conceived to overcome these limitations. 

Using the a-p scheme as an example, major differences between the current method 
and those mentioned above were explained in Sec. 2. This explicit scheme has the unusual 
property that its stability is limited only by the CFL condition, i.e., it is independent of 
p. Also, it was shown that its amplification factors are identical to those of the Leapfrog 
scheme if p = 0, and to those of the DuFort-Frankel scheme if a = 0. These coincidences 
are rather unexpected because the a-p scheme and the above classical schemes are derived 
from completely different perspectives, and the current scheme does not reduce to the 
above classical schemes in the limiting cases. 

The inviscid a-p scheme is neutrally stable and reversible in time. It is well known that 
a neutrally stable numerical analogue of Eq. (2.22) generally becomes unstable when it is 
extended to model the Euler equations. It is also obvious that a scheme that is reversible 
in time cannot model a physical problem that is irreversible in time, e.g., an inviscid flow 
problem involving shocks. Thus, the inviscid version was modified in Sec. 3 to form the a-e 
scheme. This new scheme has the unusual property that numerical diffusion is controlled 
by an adjustable parameter e. As a matter of fact, for all wavelengths, numerical diffusion 
can be uniformly bounded from above by an arbitrary small number by choosing a small 
enough e. Stability of the a-e scheme is limited by the CFL condition and 0 < e < 1. 
Moreover, if e = 0, the amplification factors of the a-e scheme are identical to those of 
the Leapfrog scheme, which has no numerical diffusion. On the other hand, if e = 1 , they 
unexpectedly become identical to each other and to the amplification factor of the highly 
diffusive Lax scheme. Note that, because the Lax scheme is very diffusive and uses a mesh 
that is staggered in time, a two-level scheme using such a mesh is often associated with 
a highly diffusive scheme. The a-e scheme, which also uses a mesh staggered in time, 
demonstrates that such a scheme could be free from numerical diffusion. 

In Sec. 4, the a-e scheme was extended to become an Euler solver. This solver has the 
unusual property that numerical diffusion at any mesh point (jf, n) can be controlled by a 
set of local parameters (e m )", m = 1,2,3. As in the a-e scheme, stability of the Euler solver 
is limited by the CFL condition and the requirement that, for all (j, n), 0 < (e m )" < 1, 
m = 1,2,3. Note that an Euler solver using a mesh staggered in time is usually highly 
diffusive for a small CFL number. It was shown in Sec. 8 that the current solver is an 
exception. It can generate highly accurate shock tube solutions with the CFL number 
ranging from 0.88 to 0.022. 

In Sec. 5, the a-p scheme was extended to become a Navier-Stokes solver. Stability 
of this explicit solver is also limited only by the CFL condition. Despite the fact that it 
does not use (i) any techniques related to the high- resolution upwind methods, and (ii) any 
ad hoc parameter, it was shown in Sec. 8 that the current solver is capable of generating 
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highly accurate shock tube solutions. Particularly, shock discontinuites can be resolved 
within one mesh interval. 

A summary of the key results of the present work has been given. Behind these results 
is a continuous effort to maintain the simplicity, generality, and accuracy of the current 
method. This effort is summarized in the following remarks: 

(a) Simplicity. The current numerical framework rests upon only two basic building 
blocks, i.e., the space-time conservation and solution elements. It uses only local dis- 
crete variables. Also, the set of discrete variables in any one of the numerical equations 
to be solved is associated with a single SE or a few immediately neighboring SE’s. 
Thus, local flexibility is preserved and one needs only to deal with a very sparse matrix. 
Moreover, flux evaluation at an interface separating two CE’s requires no interpola- 
tion or extrapolation. Nor does it require the use of an ad hoc flux model. Finally, 
partly because no characteristics-based techniques are used, a numerical scheme can 
be constructed by using only the simplest approximation techniques. 

(b) Generality. A guiding principle in the design of the current method is to limit the 
use of special assumptions or techniques that would restrict its use in more general 
situations. Thus we do not use characteristics-based techniques, and we try to avoid 
using ad hoc techniques. 

(c) Accuracy. Because (i) a physical solution of the conservation laws may involve shocks 

or high-gradient regions, and (ii) an accurate numerical simulation of such a solution 
is difficult to obtain without enforcing flux conservation, the current method requires 
that a numerical solution satisfies (i) the differential form of the conservation laws 
uniformly within an SE, and (ii) the integral form over any space-time region that is 
the union of any combination of CE’s. In addition, accuracy of the current method is 
aided by treating both (u m )" and as independent variables, instead of express- 

ing («ni)" as a finite-difference approximation involving (u m )”’s of neighboring mesh 
points. The latter approach may result in poor accuracy in a high-gradient region. 
Also, accuracy is enhanced by the fact that the flux at an interface separating two 
CE’s is evaluated without interpolation or extrapolation. Moreover, because flux con- 
servation is fundamentally a property in space-time, the current unified treatment of 
space and time may also contribute to a more accurate simulation of the conservation 
laws. 

As a result of its simplicity and generality, the current framework is also highly flexible. 
This flexibility will be demonstrated in the following discussion on how to discretize steady- 
state problems using the current method. In this discussion, we shall also address the 
important issue of boundary- condition implementation. 

As a vehicle of demonstration, we consider a dimensionless form of the 2-D steady 
incompressible Navier-Stokes equations with constant viscosity coefficient [4]. Without 
any loss of generality, we can assume that the mass density = 1. Let x and y be the first 
and second coordinates, respectively, of a 2-D Euclidean space JE2. Let u x and 112 be the 
x- and y-velocities, respectively. Let 1/3 be the static presssure. Let Rei be the Reynolds 
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number. Let 


(9.1) 

(9.2) 


£X 

fi = «i. 


rV def 

fl = «2, 


f! 


rr def / n2 , 2 dui 

U = («,) +U3 - 5 ^ ir , 

r X def 1 (dux du 2 \ 

= / s =^-rT l {-^+ aTj’ 


(9.3) 


and 


,y def , \2 I ^ &U2 

h ~ M +«3 Rsl dy ■ 


Then the Navier-Stokes equations can be expressed as 


dx + dy 


m — 1,2,3. 


(9.4) 


(9.5) 


The integral form of Eq. (9.5) in E 2 is Eq. (4.6) with h m = (/^, /m)> m = 1’^, 3, being 
the mass, i-momentum, and y-momentum flux current density vectors, respectively. Note 
that, with the understanding that the current E 2 is no longer a space-time, S(V) and ds 
have the same definitions as given in Sec. 2. 

Consider the rectangular computational domain depicted in Fig. 29(a). The upper and 
lower boundaries are in contact with stationary walls, while the left and right boundaries 
are the inlet and exit planes, respectively. The domain is first divided into rectangu- 
lar regions (see fig. 29(a)). The center of a rectangular region is denoted by its coordi- 
nates ( Xj,yk ). Each retangular region is further divided into four triangular regions (see 
Fig. 29(a)). Let the interiors of the triangular regions to the east, north, west and south of 
the center (xj, yt) be solution elements, and denoted by SE(j, fc; 1), SE(j, k: 2), SE(j, fc, 3), 
and SE(j, Jt; 4), respectively (see Fig. 29(b)-(e)). The CE’s will be defined later. 

Let q = 1,2, 3,4. For any (x,y) € SE(j,fc;<f), let u m (x,y), /^(x,y), f^(x,y), 
and h m (x,y ) be approximated by u* m (x,y ]j,k;q), /m*( x > 1 / g), and 

h* (x,y ;j,k\q), respectively. The last four functions will be defined shortly. Let (x 9 , y\) 
be a point in SE (j,k',q) (see Fig. 29(b)-(e)), and 


u* m (x,y ; j, k;q) d = (u m )’ifc + (“mil’ll 1 “ x j) + ( u my) 9 j t k(y Vk) 

+ -(Umiilj ^ 1 — x j ) 2 + 2( Urn yy) 9 j,k(y — yt) + ( u mxy)j t k( X ~ X j)(y yk )’ 


(9.6) 


where (u m ) 9 - k , ( u mx 
SE (i,fc;g). They are 
du m /dy , d 2 u m /dx 2 , 


);V (%«),V M)U> ^ ar * Giants in 

considered to be the numerical analogues of the values of u 
d 2 u m /dy 2 , and d 2 u m /dxdy at (xj, yf), respectively. 
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Similarly, for m = 1, 2, 3 and q = 1, 2, 3, 4, let 


/£•(*, I ,;>,*;«) = (fm)U + - *J) + (/i,)J,*(v - »{) 

+ - *?) 2 + 5 (/ S „) J ^( l / - S '?) 2 + (/£,,)?.*(* - ** X » - rf ). 


(9.7) 


and 


/T(*.yii.*;«) = (A)5,» + (/£«)?,*(* - *J) + - y» 

+ - X V>‘ + 5(A„)J,»(y - rf) 2 + (/£„)’,*(* - **)(» - rf)- 


(9.8) 


2 w mxx/j^ -j/ ’ 2 

Also, because h m = (/m> fm)i let 


S»(*. y (/”(*. y ; j. *; «). /«”(*> v ;i, *;*))• 


(9.9) 


The expansion coefficients in Eqs. (9.7) and (9.8) can be defined in terms of those in 
Eq. (9.6). As an example, consider (/;?**)’*. It is the numerical analogue of the value 
of d 2 f%/dx 2 at (ij, y q k ). By using Eq. (9.2), d 2 fdx 2 can be expressed in terms of t»i, 
u 3 , and their derivatives. With the understanding that the numerical analogue of any 
derivative of u m higher than second order is set to zero, one can obtain the numerical 
version of the above relation by replacing each derivative with its numerical analogue. 
Using this numerical version, (f^xxYj k can defined in terms of the expansion coeffi- 
cients in in Eq. (9.6). Note that (fm*)] and ( fmt )" were defined in a similar fashion (see 
Eqs. (4.14)-(4.17)). 

From the above discussion, one concludes that the only independent discrete variables 
needed to be solved are the expansion coefficients in Eq. (9.6). Because m = 1,2,3, there 
are 18 unknowns for each SE. 

We assume that the numerical solution satisfies Eq. (9.5) uniformly within an SE, i.e., 
for all ( x,y ) 6 SE (j,k\q), 


df„(x.y,j,k-,q) df£(x,y,j,k;q) _ 
d~x + Oy 


m — 1,2, 3. 


(9.10) 


Because Eq. (9.10) is equivalent to V • h* m — 0, Gauss’ divergence theorem implies that the 
total flux of h* m leaving the boundary of any region within an SE vanishes. 

As a result of Eqs. (9.7) and (9.8), for each m the expression on the left side of 
Eq. (9.10) is a polynomial of first order in (i — Xj) and (y — y£). Eq. (9.10) requires that all 
three coefficients of this polynomial vanish. Because m = 1,2,3, Eq. (9.10) represents nine 
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conditions for each SE. Thus to match the number of conditions with that of unknowns, 
nine more conditions per SE are needed. One of many ways to fulfill this need is described 
in the following discussion. 

First we consider an interior SE, i.e., an SE with each of its three edges being an 
interface separating two SE’s. Hereafter, an edge separating two SE’s will be referred to as 
an interior edge; while an edge bordering with the boundary of the computational domain 
will be referred to as an boundary edge. An example of an interior SE is that depicted 
in Fig. 29(f). Each interior edge is divided into two subsections. Thus (i) the edge BC 
is divided into BA' and A'C, (ii) the edge CA is divided into CB' and B'A, and (iii) the 
edge AB is divided into AC' and C'B. Moreover, for each m, we shall assume that the net 
flux of h ^ entering each subsection from both sides vanishes. Because m = 1,2,3, there 
are six interface flux conservation conditions for each interior edge. Moreover, because an 
interior edge is bordered by two SE’s, each SE can be allocated only three net interface 
flux conditions for each interior edge. Because an interior SE has three interior edges, the 
need for nine extra conditions is fulfilled by the above interface flux conditions. 

Next consider a boundary SE, i.e., an SE with at least one boundary edge. The edges 
DE and FD of the SE depicted in Fig. 29(g) axe interior edges, while the edge EF is a 
part of a stationary wadi. Again we divide the interior edges into two subsections, i.e., (i) 
DE is divided into DF' amd F'E, and (ii) FD is divided into FE' and E'D. The interface 
flux conditions imposed on DE or FD axe similar to those described earlier for other 
interior edges. As will be explained shortly, three boundary conditions will be imposed 
on a boundary edge, such ais EF. Because three net interface conditions are allocated to 
each interior edge, adding the number of boundary conditions to the number of interface 
conditions results in the nine extra conditions each boundary SE needs. This conclusion 
is valid even if the SE has more than one boundary edge. 

With the above preparations, the CE’s may be defined as follows: An interior SE such 
as that depicted in Fig. 29(f) can be divided into four triangular regions, i.e., AC'B', BA'C', 
CB'A', and A'B'C'. The union of each of these regions and its boundary, by definition, 
is a CE. On the other hand, a boundary SE such as that depicted in Fig. 29(g) can be 
divided into a triangular region DF'E' and a quadrilateral region EFE'F'. The union of 
each of these two regions and its boundary, by definition, is also a CE. 

As a result of the above definition, an interface separating two CE’s may be of two 
different types. Those of the first type Eire located within an SE. They are exemplified by 
A'B' in Fig. 29(f), and E'F' in Fig. 29(g). Because h* m is continuous on the neighborhood of 
an interface of this type, the net flux of h* m entering such an interface vanishes. Interfaces 
of the second type are subsections separating two SE’s. They are exemplified by AC' in 
Fig. 29(f), and DF' in Fig. 29(g). By assumption, the flux of entering an interface 
of this type also vanishes. Because the interior of a CE is within an SE, according to a 
statement made following Eq. (9.10), the total flux of h ^ leaving the boundary of a CE 
vanishes. Combining the results established above, one concludes that the total flux of 
leaving the union of any combination of CE’s vanishes. Note that the SE’s and CE’s 
defined here are substantially different from those depicted in Figs. 2 and 3. Furthermore, 
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it should be emphasized that it is possible to define two different sets of CE’s with the 
same SE’s and numerical conditions. 

In the above discussions, we have assumed that three boundary conditions are imposed 
on each boundary edge. Using the boundary edge EF depicted in Fig. 29(g) as an example, 
we shall describe how these boundary conditions may be implemented. 

The edge EF is a part of a stationary wall. Thus, we assume that 

ui = 0, and U 2 = 0, (9.11) 


on EF. Because EF is aligned with the r-axis, we also have 


dui 

dx 


= 0 , 


(9.12) 


on EF. Moreover, Eq. (9.12) coupled with the equation in Eq. (9.5) with m = 1 implies 
that 

-jr- = 0, (9.13) 

oy 

on EF. With the aid of Eqs. (9. 1 1)— (9.13), Eqs. (9.1), (9.2) and (9.4) imply that 


ft = f! = /a - ft = 0 , (9-14) 

on EF. As a result, we may assume that the simple average of each of ff*, ff* and 
(ff* — ff*) over the length of EF vanishes. The above assumption imposes three boundary 
conditions on EF. Note that the net numerical mass flux entering the upper wall through 
EF vanishes if the average of ff* over EF vanishes. 

At this juncture, it should be emphasized that the boundary conditions proposed 
above, as in the case of the interface flux conservation conditions, are conditions over a 
domain. For the special case under consideration, the domain is the entire length of EF. 
In general, however, the domain may be a subsection of EF. As an example, the edge EF 
may be divided into two subsections. We may require that (i) the average of ff* over 
each of these subsections vanishes, and (ii) the average of ff* over the entire length of EF 
vanishes. The resulting three conditions may replace the three conditions proposed earlier. 
Obviously, there axe many other alternatives. Futhermore, as need arises, one can easily 
impose more boundary conditions over a boundary edge. 

In the above numerical discretization, u m (x,y) is approximatd by a polynomial of 
second order within an SE (see Eq. (9.6)). Next we briefly consider other cases in which 
polynomials of first order and third order are used. 

Let the numerical version of u m (x,y) be a polynomial of first order. Then there 
are three expansion coefficients, i.e., three independent unknowns, for each m = 1,2,3. 
Thus, there are nine independent unknowns for each SE. Let the numerical solution satisfy 
Eq. (9.5) uniformly within an SE, i.e., an equation similar to Eq. (9.10) is imposed for 
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each m. Because the expression on the left side of Eq. (9.10) would become a constant 
if polynomials of first order were used in Eqs. (9.6)-(9.8), the last requirement represents 
one condition per SE for each m. Thus, to match the number of conditions with that of 
unkno wns, six extra conditions per SE are needed. 

In Fig. 30, the computational domain depicted in Fig. 29(a) is divided into triangular 
regions and rhombic regions. The interior of each triangular region is a boundary SE while 
the interior of each rhombic region is an interior SE. For each let one interface flux 
conservation condition be imposed on each interior edge. In addition, as shown earlier, 
one may impose three boundary conditions over each boundary edge. Using arguments 
given previously, it is easy to show that these interface and boundary conditions result in 
the needed six extra conditions per SE. Let the union of each SE and its boundary be a 
CE. Then one concludes that, for each m, the total flux of h* m leaving the union of any 
combination of CE’s vanishes. 

Next we consider the case in which the numerical version of y) is a polynomial of 

third order. Then there are ten expansion coefficients, i.e., ten independent unknowns, for 
each m = 1, 2, 3. Thus, there are 30 independent unknowns per SE. Again we assume that 
the numerical solution satisfies Eq. (9.5) uniformly within an SE. Because (i) the expression 
on the left side of Eq. (9.10) would become a polynomial of second order if polynomials of 
third order were used in Eqs. (9.6)-(9.8), and (ii) there are six expansion coefficients for 
each polynomial of second order involving two unknowns, the last requirement represents 
six conditions per SE for each m. Thus, to match the number of conditions with that of 
unknowns, 12 extra conditions per SE are needed. Because the number of extra conditions 
needed per SE is twice that needed in the special case we have just discussed, it is obvious 
that the current need can be fulfilled if (i) for each h* m , two interface flux conditions are 
imposed on each interior edge depicted in Fig. 30, and (ii) six boundary conditions are 
imposed on each boundary edge. To give a definition of CE’s, in Fig. 30, an interior SE 
is divided into 16 subregions, and a boundary SE is divided into six subregions. Each of 
these subregions can be considered as a CE. 

Several variants of numerical discretization for the same flow problem have been de- 
scribed. Each variant represents a system of nonlinear equations, and is implicit in nature. 
Each can be solved by a variety of solution procedures. In [3], using Newton’s method and 
another variant of discretization, a new efficient procedure was developed for the solution 
of incompressible, laminar channel flow. It was shown that, for a flow with Rti, = 100, an 
accurate solution can be obtained by using as few as six SE s across the channel. 


60 



Appendix A. An Alternative Stability Analysis for 
the Lax and Leapfrog/DuFort-Frankel Schemes 


With the use of the regular mesh depicted in Fig. 31, the Lax scheme for solving 
Eq. (2.22) can be expressed as 


~ ( u j'+i + U j'-i ) / 2 


V+i ~ ”, 

2ax' 


(A.l) 


where j',n' = 0, ±1,±2, — The system of equations represented by Eq. (A.l) can be 
divided into two sets completely independent from each other. The first set involves 
only the variables associated with those mesh points marked by dots in Fig. 31, and the 
second set, by crosses. Thus, the solution to Eq. (A.l) contains two decoupled solutions. 
Traditionally the von Neumann stability analysis for the Lax scheme is performed without 
taking into account this decoupling nature. Consider a solution to Eq. (A.l) in which 
u", = 1 for all mesh points ( j',n ') that are marked by dots, and u”, = — 1 for all other 
(j', n'). In reality, this solution represents the union of two completely decoupled constant 
solutions. However, at any time level, the combined solution is represented by a Fourier 
component of the shortest wavelength (= 2ax') in the traditional analysis. Therefore, two 
decoupled constant solutions may be wrongly perceived as a rapidly-varying solution. For 
the above reason, we shall consider each decoupled solution separately in the following von 
Neumann stability analysis. 

Let n = n'/2, j = j’ /2, ax = 2ax', and At = 2a t 1 . Then the mesh depicted in 
Fig. 31 is identical to that depicted in Fig. 2(a) except that those mesh points marked by 
crosses in Fig. 31 have no counterparts in Fig. 2(a). As a result, the decoupling nature of 
Eq. (A.l) will be removed if the Lax scheme is expressed using the staggered mesh depicted 
in Fig. 2(a), i.e., for all (j, n) € ft, 


-L ,, n_1 / 2 Wo n- 1/2 n— 1/2 

,“>+1/2 +V 1 / 2 J / 2 , _ Vn/2 -“>-1/2 

, 1 & 


With the aid of Eq. (2.13), Eq. (A.2) can be simplified as 

«? = (1/2) [(1 + ■')«;-■/’ + (1 - ^)«”- 1 1 / 2 2 ] ■ (A. 3) 

By applying Eq. (A. 3) successively, one hats 

“” +I = (1/4) [(1 + 0 2 “"-i + 2(1 - + (1 - *) 2 < + ,] . (A4) 

In contrast to Eq. (2.19), Eq. (A. 4) implies that u” +1 does not approach u” as At — ► 0. 
Moreover, by substituting 

u] = [G(u,9)] n e ij0 (i d = \/-f, -tt < 9 < tt) (A. 5) 
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into Eq. (A. 4), one concludes that the amplification factor of the Lax scheme is given by 


G(v,9) — [cos(0/2) — zVsin(0/2)] 2 . 


(A.6) 


A comparison among Eqs. (3.12), (3.15), and (A.6) reveals that — G{u,6) 

when e = 1. 

Because ti" +1 does not approach u" as At -» 0. It follows from Eq. (A.5) that G(u, 6) 
cannot approach 1 as v — ► 0. As a matter of fact, G(v,8) — > cos 2 (0/2) as u — * 0. In turn, 
this implies that the Lax scheme is highly diffusive when \v\ is small. 

With the use of the regular mesh depicted in. Fig. 31, the Leapfrog/DuFort-Frankel 
scheme for solving Eq. (2.1) can be expressed as 


u 


n'+l 


u 


t '-l 


y 


2 At' 


+ o, 


u 


y + 1 


— u 


y - 1 


2ax' 


- M 


,.n' i u n' _ -V+l 

u j>+i' u y—i *~y 




(ax') 2 


= 0. 


(A. 7) 


where j\ n' = 0, ±1, ±2, Even though Eq. (A.7) is a three-level scheme while Eq. (A.l) 

is a two-level scheme, they have the same decoupling nature. The decoupling of Eq. (A.7) 
can be removed if the scheme is expressed with respect to the staggered mesh depicted in 
Fig. 2(a), i.e., for all (j, n) €E ft, 


„ n— 1/2 n — 1/2 n-1/2 n-1/2 „ n- 1 

u j~ u j , “>+ 1/ 2 “>-1/2 “>+1/2 + U >-l/2 U i > _ n 

' +a M (ax/2) 2 


At 


AX 


With the aid of Eq. (2.13), Eq. (A.8) can be simplified as 

(1 + £)u” = (1 — £)«" 1 + {v + 0 u j- 1/2 — ( u — 0 u j+ 1/2 j 
Eq. (A. 9) can also be expressed in a two-level form, i.e., 

tf(j, n) — L+u(j — 1/2, n — 1/2) + L-u{j + 1/2, n — 1/2). 


Here 


n del 

u(;,n) = 


u. 


n — 1/2 
“>+1/2 


for all (j, n) € ft with n > 0, and 

v+t i-e 


def 


L , = 1 + i l + e , and L_ = 


0 


0 


-(*-0 

i+e 

i 


0 

0, 


(A.8) 


(A. 9) 


(A-10) 


(A. 11) 


(A. 12) 
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By applying Eq. (A. 10) successively, one has 

u(j, n + 1) = (L + ) 2 u (j - 1, n) + ( L+L _ + L-L + )u(j,n ) + (X_) 2 u(j + 1, n). (A.13) 


To perform the von Neumann stability analysis for Eq. (A. 13), let 

u(j,n) = u*(n,6)e x ^ 0 ( i *== y/— 1, — it < 6 < 7r) (A.14) 

where u*(n, 6) is a 2 x 1 column matrix. Substituting Eq. (A. 14) into Eq. (A. 13), one 
obtains 

u*(n + 1,6) = [L{ v,£,0)] 2 tf*(n, 6) (A.15) 

where 

L{v, i, 6) d = e~ i0/2 L + + e‘* /2 X_. (A.16) 

According to Eq. (A.15), [L(i/,£,0)] 2 is the amplification matrix. Substituting Eq. (A. 12) 
into Eq. (A. 16), one has 

2[£ cos(fl/2) — zVsin(fl/2)] (1-Qe~ ,g / 2 \ 

L{v,£,,6) = 1 + £ 1 + £ )• (A.17) 

\ e i0 > 2 0 / 

The amplification factors A± given in Eq. (2.21) are the eigenvalues of the amplification 
matrix [L(v, f , 9)} 2 . 
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Figure 2. — The SEs and CEs of type I. 
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(d) Three neighboring CEs. 
Figure 3. — The SEs and CEs of type II. 
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Figure 7.- The Euler solution (e = 0.5, a = 1, At = 0.004, CFL = 0.88). 
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Figure 14.- The Euler solution (e = 0.3, a = 2, At = 0.004, CFL = 0.8S). 
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Figure 25.- The Navier-Stokes solution {Rei = 8000). 
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Figure 28.- The Navier-Stokes solution (Rei = 20000). 
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(a) the computational domain. 



(b)SE(j, k; 1 ). (c)SE(j, k; 2 ). 



(d) SE (j, k; 3). (e)SE(j, k; 4). 



(f) An interior SE. (g) A boundary SE. 

Figure 29. — The SE's and CE's for a flow problem 
(first construction). 
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IT) O 


A sample program for solving the shock tube problem 


implicit real*8 (a-h, o-z) 

dimension q(3,1000), qn(3,1000), qx(3,1000), qt(3,1000), 
* s(3,1000), vxl(3), vxr(3) , xx(1000) 

it = 100 
dt = 0 . 4d-2 
dx = 0.1d-l 
ga = 1.4d0 
rhol = l.dO 
ul = O.dO 
pi = 1 . do 
rhor = 0 . 125d0 
ur = O.dO 
pr = O.ldO 
ic = 1 

hdt = dt/2.d0 
tt = hdt*df loat ( it) 
qdt = dt/4.d0 
hdx = dx/2.d0 
qdx = dx/4.d0 
dtx - dt/dx 
al = ga - l.dO 
a2 = 3. dO - ga 
a3 = a2/2 . dO 
a4 = 1 . 5d0*al 
q ( 1 , 1 ) = rhol 
q(2,l) = rhol*ul 

q(3,l) = pl/al + 0 . 5d0*rhol*ul**2 

itp = it + 1 

do 5 j = 1 , itp 

.q(l,j + l) = rhor 

q(2,j+l) = rhor*ur 

q(3,j+l) = pr/al + 0 . 5d0*rhor*ur**2 
do 5 i = 1,3 
qx(i,j) = O.dO 
continue 

open (unit=8, file= / for008 / ) 
write (8,10) tt,it,ic 
write (8,20) dt,dx,ga 
write (8,30) rhol,ul,pl 
write (8,40) rhor,ur,pr 

m = 2 

do 400 i = l,it 
do 100 j = l,m 
w2 = q(2, j)/q(l, j) 
w 3 = q(3, j)/q(l, j) 


C’l 
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100 


200 


300 

400 

c 


500 


600 

C 

10 

20 

30 

40 

50 


f21 = -a3*w2**2 
f 22 = a2*w2 

f31 = al*w2**3 - ga*w2*w3 
f32 = ga*w3 - a4*w2**2 
f 33 = ga*w2 
qt(l#j) = -qx(2 f j) 

qt ( 2 , j ) = - ( f 21*qx (1 , j ) + f22*qx(2,;j) + al*qx(3,3)) 
qt ( 3 , j ) = - ( f 31*qx ( 1 / j ) + f32*qx(2,j) + f33*qx(3,j)) 
s(l,j) = qdx*qx ( 1 , j ) + dtx*(q(2,j) + qdt*qt(2,j)) 
s (2 , j ) = qdx*qx(2,j) + dtx* (f 21* (q(l, j ) + qdt*qt ( 1 , d ) ) + 

* f22*(q(2,j) + qdt*qt(2,j)) + al*(q(3,j) + qdt*qt (3 , j ) ) ) 
s (3 , j ) = qdx*qx(3,.j) + dtx* (f 3 1* (q(l, j ) + qdt*qt(l,j)) + 

* f32*(q(2,j) + qdt*qt (2 , j ) ) + f33*(q(3,j) + qdt*qt (3 , j ) ) ) 
continue 

mm = m - 1 
do 200 j = l,mm 
do 200 k =1,3 

qn (k, j+1) = 0.5d0*(q(k, j) + q(k,j+l) + s(k,j) - s(k,j+i)) 

vxl(k) = (qn(k,j+l) - q(k,j) - hdt*qt(k, j) )/hdx 

vxr(k) = (q(k,j+l) + hdt*qt (k, j+1) - qn(k, j+1) )/hdx 

qx (k, j+1) = (vxl(k) * (dabs (vxr(k) ) )**ic + vxr(k)* (dabs (vxl(k) ) ) 

* **ic)/ ( (dabs (vxl (k) ) ) **ic + (dabs (vxr (k) ) ) **ic + l.d-60) 
continue 

do 300 j = 2,m 
do 300 k = 1,3 
q(k,j) = qn ( k , j ) 
continue 
m = m + 1 
continue 

t2 = dx*df loat (itp) 
xx ( 1 ) = -0 . 5d0*t2 
do 500 j = 1, itp 
xx (j+1) = xx(j) + dx . 
continue 
do 600 j = l,m 
x = q(2, j)/q(l, j) 

z = al* (q(3 , j ) - 0 . 5d0*x**2*q( 1, j ) ) 
write (8,50) xx ( j ) , q(l , j ) , x, z 
continue 


close (unit=8) 

format ( ' t = ',gl4.7,' it = ',i4, # ic = # ,i4) 

format ( ' dt = ',gl4.7, 7 dx = / ,gl4.7, / gamma = / ,gl4.7) 

format (' rhol = , ,gl4.7, / ul = / ,gl4.7, / pi = ',gl4.7) 

format ( 7 rhor = / ,gl4.7, / ur = # # gl4. 7, 7 pr = ',gl4.7) 

format ( ' x = / ,f8.4, / rho =',gl4.7, # u ="^14. 7, 7 p = / ,gl4.7) 

stop 

end 
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